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1  EXECUTIVE  SUMMARY 


1.1  Objectives 

The  objectives  of  this  program  are  aimed  at  the  development  of  controls  technology  for  Smart 
Structures  and  Materials  systems  with  large  numbers  of  sensors  and  actuators.  This  develop¬ 
ment  focuses  on  a  methodology  for  designing  controllers  that  control  large-scale  Smart  Struc¬ 
tures  and  Materials  systems  with  potentially  thousands  of  inputs  and  outputs  at  a  dramatically 
reduced  cost  in  both  computation  and  input/output  interconnection  complexity.  Moreover, 
the  methodology  can  be  used  to  implement  large-scale  controllers  that  have  highly  scalable, 
multiprocessor  architectures. 

Our  program  objectives  can  be  summarized  as  follows: 


•  Demonstrate  an  initial  proof-of-concept  of  the  wavelet-based  method  for  controlling 
smart  structures  and  materials  with  large  transducer  arrays. 

•  Develop  a  practical  and  mathematically  sound  methodology  for  designing  and  imple¬ 
menting  computationally  efficient,  sparsely  interconnected,  scalable  control  systems. 

•  Validate  the  new  methodology  by  designing  and  simulating  a  control  system  for  control 
of  simple  structures  using  both  numerical  modeling  data  and  experimentally  measured 
data. 


1.2  Results 

The  Control  Technology  for  Smart  Materials  Program  has  confirmed  the  advantages  of  a  mul¬ 
tiscale  approach  to  transfer  function-based  control  of  systems  with  very  large  numbers  of 
transducers.  The  main  results  of  this  program  are  as  follows: 


•  Wavelet  Control  Design  Method 

Main  result: 

Developed  design  method  that  exploits  wavelet  transform  compression  properties  and 
the  Q-parameterization  of  the  controller  in  producing  controller  designs  that  can  be 
implemented  with  reduced  computational  processing  demands  and  with  a  simplification 
of  interconnections  from  inputs  to  outputs. 

Summary  of  method: 

By  using  the  Q-parameterization  of  the  controller,  we  can  divide  the  design  problem  into 
the  two  parts.  The  first  part  consists  of  implementing  a  transfer  function  that  cancels  the 
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coupling  from  actuator  inputs  to  sensor  outputs.  The  second  part  consists  of  design¬ 
ing  a  feedforward  transfer  function,  “Q”,  such  that  its  parameters  minimize  an  optimal 
control  cost  function.  The  wavelet  control  method  uses  the  wavelet  transform  to  com¬ 
press  the  coupling  transfer  function,  and  uses  a  spatially  banded  Q-parameter  structure 
to  compress  the  feedforward  transfer  function.  The  result  is  an  overall  controller  that  is 
sparsely  interconnected  and  that  has  fewer  effective  input/output  interconnections. 

Our  design  method  also  exploits  temporal  bandwidth  compression  via  finite-impulse 
response  (FIR)  parameterization  of  the  controller  and  multirate  filtering. 

•  Sensor/Actuator  Placement  Methods  and  Integration  with  KIKO 

Main  results: 

-  Developed  algorithm  for  transducer  selection.  This  approach  uses  a  predefined  set  of 

candidate  actuators  or  sensors,  and  selects  a  fixed,  smaller,  number  of  these  trans¬ 
ducers  to  maximize  their  effectiveness.  For  this  purpose  effectiveness  is  defined 
for  actuators  as  approximating  the  range  of  all  possible  uncontrolled  system  out¬ 
puts  over  a  number  of  frequencies,  and  for  sensors  as  approximating  the  span  of 
all  possible  system  disturbance  inputs  over  a  number  of  frequencies. 

-  Developed  approaches  for  using  nonuniform  and  irregularly  bounded  distributions 

of  transducers.  The  result  of  transducer  selection  would  generally  be  irregularly 
placed  devices.  There  is  no  theoretical  basis  for  benefits  of  direct  application 
wavelet  transformation  in  this  case,  so  we  looked  at  three  ideas  to  retain  the  ad¬ 
vantage  of  multiscale  methods  with  irregular  spatial  placement  of  transducers. 

•  Validation  Using  Analytical  and  Numerical  models 
Main  results: 

-  Validation  of  wavelet  compression  method  on  Euler-Bemoulli  beam  model: 

Used  transfer  function  from  force  inputs  to  displacement  outputs  along  a  beam 
as  example  of  coupling  transfer  function.  Showed  compression  performance  that 
agreed  with  theoretical  predictions  over  wide  range  of  frequencies 

-  Validation  of  overall  control  method  on  Euler-Bemoulli  beam  model: 

Successfully  controlled  displacement  along  a  beam  due  to  wide-band  force  dis¬ 
turbances.  Wide-band  control  performance  achieved  using  a  wavelet-based  con¬ 
troller. 

•  Validation  Using  PSU-collected  Data 

Main  results: 

-  Validation  of  wavelet  compression  method  on  measured  plate  transfer  function: 
Used  34  X  34  X  .125  inch  aluminum  plate  to  measure  coupling  transfer  function. 

-  Validated  smoothness  properties  of  plate  using  taking  numerical  approximation  of 
derivatives  using  measured  data 

-  Validated  wavelet  compression  of  transfer  function  using  measured  data 
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1.3  Recommendations 


The  next  logical  steps  of  development  are  closed-loop  hardware  demonstration,  and  further 
work  on  advanced  wavelets  that  can  treat  irregular  transducer  arrays,  accomodate  bounded 
domains,  and  reduce  the  break-even  point  for  advantageous  application  of  this  technology. 
The  first  step  is  needed  to  build  confidence  of  potential  users,  and  is  justified  by  the  results  to 
date.  The  second  step  would  be  focussed  on  the  most  important  needs  for  further  theoretical 
and  algorithmic  development  of  the  methodology  to  suit  real  applications. 
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2  INTRODUCTION 


Smart  Structures  and  Materials  technology  has  emerged  as  a  promising  new  area  of  engi¬ 
neering  research  and  development.  The  concept  of  embedding  sensing  and  actuating  devices 
within  a  material  and/or  structure  offers  the  possibility  of  altering  the  physical  behavior  of  a 
system  in  ways  not  possible  through  the  use  of  passive  materials  alone.  The  DARPA  Smart 
Structures  and  Materials  program  has  focused  much  of  its  effort  on  the  development  of  de¬ 
vice  and  embedding  technology  for  initial  proof  of  concept  for  engineering  such  systems. 
The  design  requirements  of  embedding  devices  in  materials  and  structures  present  formidable 
challenges  to  developing  device  technology  capable  of  performing  in  applications  requiring, 
for  example,  the  rendering  of  significant  actuation  force  over  a  reasonable  bandwidth. 

The  development  of  active  control  technology  specifically  for  smart  structures  and  materials, 
however,  has  lagged  substantially  behind  that  of  the  base  materials,  transducers,  and  embed¬ 
ding  techniques.  Still,  development  of  smart  structures  with  ever-greater  numbers  of  embed¬ 
ded  elements  continues,  spurred  by  potential  uses  that  require  large  arrays  of  sensors  and 
actuators.  For  example,  several  “Smart  Skin”  concepts  are  being  developed  for  embedding 
transducers  on  the  surface  of  a  structure.  By  placing  them  on  the  wing  of  an  aircraft,  one 
increases  flight  performance  by  detailed  control  of  the  wing’s  shape,  i.e.,  a  “Smart  Wing.” 
By  embedding  transducers  in  an  aircraft  cabin’s  walls,  one  achieves  interior  acoustic  noise 
quieting.  There  is  yet  another  reason  for  considering  an  array  of  transducers  on  the  surface 
of  a  structure.  Rather  than  implementing  a  conventional  control  design  that  is  sensitive  to  the 
particular  device  layout,  a  densely  sampled  array  allows  the  controller  to  optimize  the  use  of 
sensor  information  and  actuator  authority.  No  control  technology  suitable  for  such  large  arrays 
exists,  however,  and  this  presents  a  barrier  to  future  applications. 

To  motivate  the  objectives  of  our  program,  consider  the  example  of  controlling  a  large  array 
of  collocated  actuator/sensor  pairs  on  a  flexible  plate  as  illustrated  in  Figure  1.  Each  tile 
on  the  surface  contains  an  embedded  actuator/sensor  pair  and  the  control  problem  might  be, 
for  example,  to  stiffen  the  plate  by  dynamically  adjusting  the  actuator  input  signals  based  on 
dynamically  sensing  the  sensor  output  signals.  In  Figure  1,  the  abundance  of  communication 
paths  illustrates  the  complexity  involved  in  mapping  the  sensor  outputs  on  the  surface  into  a 
control  signal  for  a  single  actuator  using  a  fully  interconnected  control  approach.  Note  that 
the  complexity  of  this  approach  in  terms  of  both  computational  and  communications  cost 
increases  with  the  number  of  actuators  being  controlled.  In  the  case  of  full  interconnection 
between  every  actuator  and  sensor,  the  system  becomes  potentially  intractable  and  highly  cost- 
ineffective. 

In  this  program  we  develop  and  demonstrate  technology  capable  of  controlling  hundreds  or 
thousands  of  sensors  and  actuators  embedded  in  the  base  material.  We  have  dubbed  this  the 
“KIKO  control  problem”  (Kilo-Input/Kilo-Output)  for  Smart  Structures  and  Materials.  This 
program  focuses  on  a  new  multiscale/multirate  theory  of  hierarchical  design  based  on  the 
wavelet  transform.  In  the  context  of  this  theory,  we  develop  efficient  and  highly  scalable 
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Figure  1:  Fully  interconnected  control  approach 

implementations  of  control  systems  using  multiprocessor  architecUires.  This  report  describes 
our  multiscale  control  approach  and  presents  simulation  results  based  on  model  and  measured 
data  that  provide  a  first-order  proof  of  concept  for  our  approach. 

Our  program  is  aimed  at  developing  a  theory  and  methodology  for  controlling  large  num¬ 
bers  of  devices  for  applications  such  as  vibration  control  and  acoustical  noise  control.  The 
methodology  we  develop  will  allow  us  to  implement  controllers  that  can  efficiently  process 
large  numbers  of  inputs  and  outputs.  It  is  highly  scalable  and  will  guide  us  in  designing  highly 
structured  multiprocessor  schemes.  The  end  result  is  a  control  system  that  implements  large- 
scale  control  at  a  highly  reduced  cost  in  terms  of  processing  hardware  and  interconnection 
complexity  compared  with  a  conventional  approach.  Our  approach  to  control  design  is  based 
on  a  multiscale  hierarchy  whose  structure  is  governed  by  the  characteristics  of  the  control 
system  transfer  functions  as  illustrated  in  Figure  2.  Aggregating  actuator  inputs  and  sensor 
outputs  into  local  groups  and  doing  this  recursively  in  scale  results  in  a  processing  architec¬ 
ture  that  allows  efficient  computation  and  communication  between  all  sensors  and  actuators. 
Note  that  in  contrast  to  locally  interconnected  control  approaches,  where  sensors  and  actuators 
are  constrained  to  communicate  within  local  neighborhoods,  our  approach  allows  direct  com¬ 
munication  from  any  sensor  to  any  actuator.  Thus,  while  the  controller  has  the  same  number 
of  degrees  of  freedom  as  a  fully  interconnected  controller,  it  exploits  the  multiscale  hierarchy 
in  significantly  reducing  computational  and  communications  cost. 

This  spatial  hierarchy  also  allows  us  to  use  a  multirate  approach  to  efficiently  implement  the 
controller.  Thus,  we  optimize  not  only  spatial  properties  of  the  system  to  reduce  intercon¬ 
nection  complexity,  but  temporal  bandwidth  as  well.  One  of  the  ultimate  features  of  such  a 
scheme  is  to  allow  the  control  of  broadband  phenomena  over  a  wide  spatial  extent  and  to  do 
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Figure  2:  Multiscale  control  approach 
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so  without  having  to  optimize  the  layout  of  a  relatively  small  number  of  transducers. 

This  report  sununarizes  the  theory  developed  under  this  program  and  describes  the  multiscale 
method  for  controlling  smart  structures  and  materials  in  Sections  4  and  5.  Section  1  provides 
an  executive  summary  of  our  results.  Sections  6  and  7  provide  examples  of  simulations  using 
analytical  and  numerical  models  as  well  as  results  of  applying  our  method  to  measured  data. 
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3  CONTROL  PROBLEM  FORMULATION 


This  section  describes  the  general  closed-loop  control  problem  in  terms  of  a  block  diagram  and 
its  associated  transfer  function  equations.  The  problem  formulation  applies  to  a  wide  range  of 
control  problems  and  provides  us  with  basic  notation  and  a  formalism  for  analyzing  closed- 
loop  control  systems  in  terms  of  stability  and  robustness.  Furthermore,  by  parameterizing  the 
controller  via  the  Q-parameterization  [9],  the  controller  design  problem  becomes  a  tractable 
model-matching  problem.  SRI  is  a  pioneer  in  exploiting  the  Q-parameterization  for  the  design 
of  large-scale  systems.  The  KIKO  control  approach  extends  this  approach  much  further  by 
exploiting  inherent  redundancies  in  the  physical  system  to  achieve  dramatic  reductions  in  the 
computational  and  communications  complexity  of  the  controller. 

The  block  diagram  for  the  basic  closed-loop  control  configuration  is  depicted  in  Figure  3. 

The  quadruple  of  transfer  functions  (T,  S,  A,  P)  represents  various  sub-blocks  of  the  overall 
physical  system.  In  the  absence  of  control,  the  disturbance  propagates  through  the  system  T  to 
produce  a  response,  while  the  transfer  function  K  represents  the  user-designed  controller.  The 
transfer  function  T  represents  the  uncontrolled  response  that  is  to  be  “matched”  or  canceled 
by  the  control  path.  The  transfer  function  S  maps  the  disturbance  to  the  output  of  the  sensors, 
while  A  maps  the  inputs  to  the  actuators  to  the  physical  response  variable  to  be  controlled. 
Finally,  P  represents  the  coupling  from  actuators  to  sensors.  The  overall  closed-loop  transfer 
function  in  Figure  3  from  input  to  output,  which  we  denote  G,  is  expressed  as  follows. 

G  =  T  +  A{I-KP)-^KS  (1) 

where  each  transfer  function  is  a  matrix  function  of  frequency  and  /  denotes  the  identity 
matrix,  which  is  constant  for  all  frequencies. 

Control  problems  can  be  stated  in  terms  of  an  optimization  problem  in  which  the  optimal 
controller  is  the  solution  to  the  minimization  of  the  norm  of  the  transfer  function  G  with 
respect  to  K. 

imn||G|l  (2) 

Constraints  to  this  optimization  problem  include  requiring  stability  of  G  as  well  as  insuring 
robustness  of  the  design  with  respect  to  modeling  errors  in  (T,  S,  A,  P).  Furthermore,  K  must 
be  causal,  i.e.,  realizable.  An  obvious  difficulty  in  solving  such  a  constrained  optimization  is 
the  nonlinear  dependence  of  G  on  K. 

Fortunately,  we  can  simplify  the  overall  optimization  problem  with  the  use  of  the  Q-parameterization 
of  the  controller  K.  Consider  the  set  of  controllers  given  by  the  formula 

K  =  (I  +  QP)-'Q  (3) 

where  K  is  a.  free  transfer  function  parameter.  Such  a  system  K  has  a  natural  implementation 
as  a  multi-input/multi-output  (MIMO)  feedback  system  as  depicted  in  Figure  4.  The  basic 
theorem  resulting  from  this  parameterization  is  the  following  [9]: 
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Figure  3:  Block  diagram  of  basic  closed-loop  system 


Figure  4:  Block  diagram  of  Q-parameterized  controller 


•  Referring  to  Figure  3  and  eq.(3),  for  stable  P,  the  collection  of  stable  Q  ’s  parameterize 
the  collection  of  K’s  that  stabilize  the  closed-loop  system  G. 


The  design  implications  of  this  theorem  are  far-reaching.  In  particular,  by  substituting  eq.(3) 
into  eq.(l),  we  get  the  following  expression  for  G: 

G  =  r  -t-  AQS  (4) 

Essentially,  we  have  transformed  G  from  a  nonlinear  function  of  K  into  an  affine  function  of 
Q.  Furthermore,  the  overall  optimization  problem  as  stated  in  eq.(2)  reduces  to  the  following 
model  matching  problem  [9]. 

mrailr  +  ^aSII  (5) 

If  we  restrict  Q  to  be  stable,  then,  via  the  Q-parameterization  theorem,  we  satisfy  the  require¬ 
ment  that  the  overall  closed-loop  system  G  be  stable. 

For  practical  reasons,  we  would  like  to  design  our  controller  to  be  robust  to  unmodeled  errors 
in  the  overall  system.  Indeed,  one  of  the  main  reasons  for  designing  feedback  systems  is  to 
make  a  system  insensitive  to  modeling  errors  or  changes  in  plant  parameters.  For  example, 
consider  the  case  where  we  have  additive  plant  errors.  In 

P  =  P  ad  (6) 

P  represents  the  nominal  plant  and  D  is  a  stable  transfer  function  that  represents  a  frequency- 
dependent  additive  error.  In  this  case,  the  loop  gain  of  our  system,  G,  is  equal  to  DQ,  and  by 
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the  small  gain  theorem  [8],  the  following  is  a  sufficient  condition  for  stability  of  G  given  the 
nominal  plant  P  and  the  plant  errors  as  in  eq.(6). 

\\DQ\\oo  <  1  (7) 

where  the  oo-norm  of  a  stable  matrix  transfer  function  is  the  supremum  of  the  maximum  sin¬ 
gular  value  of  the  matrix  over  all  frequencies.  A  version  of  the  robust  control  design  problem, 
then,  is  using  the  nominal  plant  P  as  our  plant  model  and  solving  eq.(5)  for  Q  under  the 
constraint  given  in  eq.(7)  as  D  varies. 

Finally,  we  consider  the  computational  complexity  of  implementing  a  controller  as  given  in 
eq.(3).  One  of  the  main  difficulties  in  implementing  a  controller  with  large  numbers  of  sensors 
and  actuators  is  the  fact  that  we  must  map  dynamically  all  of  the  sensor  outputs  into  all  of  the 
actuator  inputs.  In  the  Q-parameterization  realization  of  the  controller,  the  implementation 
problem  can  be  divided  into  two  steps: 

Step  1.  Implementing  a  model  of  P  as  a  transfer  function  from  actuator  inputs  to  sensor 
outputs  (sometimes  referred  to  as  the  neutralization  path). 

Step  2.  Implementing  (5  as  a  transfer  function  from  sensor  outputs  to  actuator  inputs  (some¬ 
times  referred  to  as  the  feedforward  control  path). 

Thus,  if  P  and  Q  are  full  matrices  over  a  wide  band  of  frequencies,  the  computational  com¬ 
plexity  of  implementing  a  controller  with  large  numbers  of  sensors  and  actuators  becomes 
significant.  For  example,  suppose  P  and  Q  were  implemented  as  matrix  finite  impulse  re¬ 
sponse  (FIR)  filters  and  there  were  Ns  sensors  and  Na  actuators.  The  number  of  multiplica¬ 
tions  and  additions  required  to  perform  a  single  tap  filter  operation  is  2NsNa.  Moreover,  the 
communications  complexity  scales  poorly,  as  illustrated  in  Figure  1. 

Finally,  the  computational  complexity  scales  linearly  with  the  number  of  taps  of  each  FIR 
filter,  motivating  the  need  for  multirate  processing.  Our  multiscale  hierarchy  gives  rise  to  a 
controller  in  which  some  channels  are  sampled  at  a  lower  rate,  while  others  are  processed 
using  shorter  FIR  filters. 
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4  WAVELET  CONTROL  DESIGN  METHOD 


4.1  Overview 

In  this  section  we  describe  in  detail  the  JQKO  control  design  method.  The  method  is  based 
on  the  Q-parameterization  of  the  controller  and  the  wavelet  decomposition  of  particular  sys¬ 
tem  transfer  functions  in  this  parameterization.  Referring  to  the  block  diagram  in  Figure  4, 
the  controller  can  be  designed  and  implemented  in  terms  of  the  feedforward  transfer  function, 
“Q,”  and  the  feedback  or  neutralization  path,  “R”  Motivated  by  our  twin  objectives  of  re¬ 
ducing  computational  and  communications  complexity  of  the  controller,  we  develop  methods 
for  designing  the  two  transfer  function  paths  such  that  they  are  both  sparse  and  have  regu¬ 
lar  structure.  Sparsity  is  a  mathematical  property  that  translates  into  reduced  computation  in 
implementing  the  controller,  while  regularity  of  the  input/output  mapping  of  the  transfer  func¬ 
tions  can  be  exploited  to  build  multiprocessor  architectures  that  are  interconnected  in  a  simple, 
scalable  fashion. 


4.2  Change  of  Basis  for  Control  Problem 

In  this  section  we  address  the  problem  of  simplifying  the  computations  required  to  implement 
a  controller  K  via  a  change  of  basis  in  both  the  space  of  sensor  outputs  (i.e.,  the  domain  of 
the  transfer  function  S)  as  well  as  the  space  of  actuator  inputs  (i.e.,  the  domain  of  the  transfer 
function  A).  The  main  idea  behind  reducing  the  computational  requirements  of  the  controller 
via  a  change  of  basis  consists  of  (1)  using  basis  transformations  to  decouple  both  the  neutral¬ 
ization  and  feedforward  transfer  functions,  i.e.,  making  the  transformed  transfer  functions  P 
and  Q  sparse,  and  (2)  using  basis  transformations  that  are  computationally  efficient  and  struc¬ 
tured  in  such  a  way  as  to  be  easily  implementable  in  digital  hardware.  In  the  following  we 
define  notation  and  a  formalism  for  describing  and  analyzing  such  an  idea. 

From  the  previous  section,  and  in  particular  from  the  parameterization  of  the  controller  K  as 
given  in  eq.(3),  there  are  two  matrix  transfer  functions  for  which  we  would  like  to  consider 
basis  transformations  as  a  way  of  simplifying  implementation,  namely,  P  and  Q.  For  our 
control  applications,  P  and  Q  are  in  general  full  transfer  function  matrices  over  the  band  of 
interest,  necessitating  fully  interconnected  systems  for  carrying  out  steps  1  and  2  from  the 
previous  section.  The  intent  of  performing  basis  transformations  on  these  transfer  functions  is 
to  sparsify  the  required  interconnections  via  highly  structured,  computationally  efficient  basis 
transformations.  We  restrict  ourselves  to  the  case  of  frequency-independent  bases  (constant 
matrices).  In  particular,  we  define  the  following: 


p  =  XiPX;^ 

(8) 

Q  =  YiQY-^ 

(9) 
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Note  that  for  both  P  and  Q,  separate  basis  transformations  are  performed  for  the  domain 
and  range  of  each  respective  matrix  transfer  function.  Indeed,  neither  P  nor  Q  is  necessarily 
square.  Note,  however,  that  in  the  case  of  square  P  {Q),  if  X~^  =  X~^  {Y~^  —  Yp^),  then  X 
and  Y  can  be  viewed  as  similarity  transformations. 

The  efficiency  of  a  change  of  basis  in  reducing  computational  complexity  depends  on  (1)  the 
sparseness  of  P  and  Q,  i.e.,  the  interconnections  required  in  realizing  P  and  Q,  and  (2)  the 
computational  complexity  of  the  basis  change  transformation.  To  see  this  more  explicitly, 
consider  that  the  implementational  requirements  of  a  system  implemented  in  bases  X  and  Y 
can  be  divided  into  the  following  two  tasks,  which  correspond  to  the  two  tasks  as  defined  in 
the  previous  section  for  implementing  the  controller  K. 


Neutralization: 

Xf^PXr 

(10) 

Feedforward: 

Yr^QY 

(11) 

In  general  it  is  not  clear  to  what  extent  one  can  sparsify  P  and  Q,  i.e.,  render  basis  coefficients 
precisely  equal  to  zero,  by  using  basis  transformations  X  and  Y  that  are  computationally  ef¬ 
ficient  to  implement.  For  example,  Jordan  and  Schur  canonical  forms  [10]  offer  similarity 
transformations  of  a  matrix  into  nearly  diagonal  and  triangular  forms,  respectively.  In  both 
cases,  however,  full  transformation  matrices  are  involved  with  0{N^)  computational  com¬ 
plexity  required  for  performing  the  transformations. 

Our  approach  to  sparsifying  the  control  transfer  functions  P  and  Q  using  basis  transformations 
is  to  develop  a  methodology  for  trading  off  sparsity  of  the  transformed  transfer  functions  with 
both  computational  and  implementational  complexity  of  the  transformations.  To  develop  this 
methodology  we  need  tools  for  characterizing  this  tradeoff  in  terms  of  analytical  properties 
of  P  and  Q.  In  particular,  we  base  our  methodology  on  tools  developed  in  the  context  of 
wavelet  transform  theory  [17].  Within  the  context  of  wavelets,  as  will  be  described  in  the 
following  section,  we  have  available  a  class  of  extremely  fast  (0{N))  basis  transforms  that 
have  the  added  capability  of  sparsifying  physical  operators.  Furthermore,  the  degree  to  which 
an  operator  can  be  made  sparse  can  be  characterized  in  terms  of  mathematical  characteristics 
of  the  operator,  such  as  its  smoothness. 

In  summary,  our  approach  to  multiscale  control  is  to  perform  multiscale  basis  transformations 
on  P  and  Q  and  implement  the  controller  in  the  space  of  multiscale  basis  coefficients.  We 
design  our  transforms  so  as  to  make  the  basis  coefficients  as  decoupled  as  possible.  This  will 
necessarily  involve,  as  will  be  explained  in  the  following  sections,  approximations  of  P  and 
Q.  However,  we  will  be  able  to  arbitrarily  trade  off  approximation  precision  with  sparsity  in  a 
structured  way. 
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4.3  Miiltiscale  Bases 


In  this  section  we  describe  an  approach  to  multiscale  control  that  is  based  on  a  class  of  ba¬ 
sis  transformations  corresponding  to  the  compactly  supported,  orthonormal  wavelet  transform 
described  in  Daubechies  [7].  The  advantages  of  using  this  particular  class  of  transforms  for 
performing  basis  transformations  on  the  plant  P  and  the  control  parameter  Q  are  threefold: 
(1)  these  transforms  take  advantage  of  scale-recursive  algorithms  that  are  extremely  efficient 
{P{N),  where  N  is  the  number  of  samples  of  the  discrete  signal  to  be  transformed);  (2)  the 
basis  functions  can  be  designed  to  be  extremely  effective  in  compressing  operators,  e.g.,  re¬ 
ducing  the  interconnectivity  of  transfer  functions  such  as  P  and  Q\  and  (3)  by  doing  a  spatial 
scale  decomposition  using  wavelets,  we  actually  perform  an  approximate  modal  decomposi¬ 
tion  that  allows  us  to  process  fine-scale  coefficients  (representing  higher  modes  that  are  more 
highly  damped)  using  shorter  impulse  responses,  and  coarse-scale  coefficients  (representing 
lower  modes  that  are  of  narrower  bandwidth)  using  lower  sampling  rates.  Section  4.3.1  gives 
a  brief  overview  of  orthonormal  wavelet  bases.  Section  4.3.2  discusses  results  of  using  these 
bases  to  sparsify  operators. 


4,3.1  Orthonormal  Wavelets  of  Compact  Support 

The  basic  structure  of  any  scale-based  wavelet  basis  is  based  on  the  property  that  a  space 
of  functions,  for  example  can  be  spanned  by  a  basis  set  consisting  of  dilations  and 

translations  of  a  single,  properly  chosen  function. 

{^(2*” a:  —  n)  for  m,  n  integer}  (12) 

The  basic  wavelet  function  in  eq.(12)  is  in  general  oscillatory  and  compressed  in  x,  hence 
the  name.  For  our  purposes  one  of  the  most  interesting  aspects  of  bases  of  the  type  in  eq.(12) 
is  the  fact  that  under  certain  conditions,  transforming  a  signal  in  these  basis  functions  can 
be  interpreted  as  performing  what  is  known  as  a  multiresolution  analysis  of  the  signal  [12]. 
In  this  analysis,  the  signal  is  represented  at  different  resolutions  using  a  scaling  function,  $, 
which  is  used  to  blur  the  signal  at  at  these  resolutions. 

The  multiresolution  analysis  interpretation  of  the  wavelet  transform  can  be  summarized  by  the 
following  key  points. 


•  The  translations  of  the  scaling  function  for  a  fixed-scale  m,  {$(2’"a;  —  n)  for  n  integer}, 

define  a  subspace  of  Representing  a  signal  at  this  scale  is  achieved  by  projecting 

the  signal  onto  this  subspace. 

•  The  translations  of  the  wavelet  function  for  a  fixed-scale  m,  {^(2™x-n)  for  n  integer}, 
span  the  subspace  defined  by  the  difference  between  the  scaling  subspace  at  scale  m  and 
the  scaling  subspace  at  scale  m  —  1. 
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c{m-2,n) ,  d(m-2,n) 


c(m-1,n) ,  d(m-1,n) 
c(m,n) ,  d(m,n) 


Figure  5:  Domain  of  wavelet  coefficients  on  a  multiscale  lattice 


•  The  wavelet  basis  coefficients  can  be  computed  recursively  in  scale,  resulting  in  an 
extremely  efficient  procedure. 


In  Daubechies[7],  extensions  of  the  Haar  wavelet  are  derived  for  a  class  of  wavelets  that  are 
compactly  supported  and  orthonormal.  Furthermore,  a  characteristic  of  this  class  of  wavelets 
is  that  their  transforms  can  be  computed  efficiently  via  a  scale  recursion  as  illustrated  in  the 
Haar  case.  Let  us  define  the  scaling  coefficients  and  wavelet  coefficients  as  follows  for  some 
signal  f{x). 

/OO 

-OO 

/OO 

-OO 

It  can  be  shown  [7]  that  the  following  recursion  can  be  used  to  compute  the  scaling  and  wavelet 
coefficients  recursively  in  scale  from  fine  to  coarse. 


c(m,  n) 

=  -{■\,k)h{2n  —  k) 

k 

(15) 

d{rn,n) 

=  Yh  1  >  ^)5'(2n  -  k) 

(16) 

k 


Each  of  the  recursions  in  eqs.(15,16)  can  be  interpreted  as  discrete  filtering  followed  by  down- 
sampling  by  a  factor  of  two.  The  scaling  function  filter  h{n)  is  related  to  the  wavelet  filter  g{n) 
by:  g{n)  =  (— 1)”/j(1  -  n).  Furthermore,  for  compactly  supported  wavelets,  h{n)  and  g{n) 
are  FIR  filters  (typically  of  length  <  12).  The  computational  complexity  for  the  algorithm  de¬ 
scribed  by  eqs.(15,16)  is  order  0{NL),  where  L  is  the  length  of  h  and  g,  and  N  is  the  number 
of  coefficients  at  the  initial  (finest)  scale.  The  algorithm  flow  is  illustrated  in  Figure  5  for  the 
case  of  4-tap  FIR  filters  h  and  g,  where  the  c  and  d  coefficients  are  computed  recursively  up 
the  lattice.  The  algorithm  is  initialized  at  the  finest  scale  with  the  sequence  c(M,  n),  where  M 
is  the  finest  scale  and  c{M,n)  is  the  signal  to  be  transformed. 

The  inverse  transform,  i.e.,  reconstructing  the  scaling  coefficients  at  the  finest  scale  from 
wavelet  coefficients  at  all  coarser  scales,  consists  of  the  following  recursion  in  scale  from 
coarse  to  fine. 

c(m  +  1,  n)  =  ^  c{m,  k)h{n  —  2k)  +  ^  c?(m,  k)g{n  —  2k)  (17) 

k  k 


)f{x)dx 

(13) 

)f{x)dx 

(14) 
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sensor  outputs 

Figure  6:  Multiscale  basis  transformation  for  control 

The  algorithm  flow  would  correspond  to  propagating  down  the  lattice  in  Figure  5.  Again,  the 
computational  complexity  is  0{NL). 

What  these  fast  wavelet  algorithms  suggest  is  the  possibility  for  performing  basis  transfor¬ 
mations  on  our  control  transfer  functions  P  and  Q,  and,  if  the  basis  coefficients  behave  in  a 
decoupled  fashion,  for  performing  local  computations  on  each  basis  coefficient  to  implement 
the  control  law.  Furthermore,  the  lattice  algorithm  suggests  a  scheme  for  parallel  processing  in 
which  processors  are  assigned  to  each  node  in  the  lattice.  Figure  6  illustrates  a  block  diagram 
of  the  flow  of  information  involved  in  implementing  the  feedforward  path  of  the  controller. 
The  neutralization  path  is  similar  but  runs  in  the  opposite  direction,  i.e.,  from  actuator  inputs 
to  sensor  outputs. 
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4.3.2  Compression  Properties  of  Waveiets 


In  this  section  we  review  some  results  in  Beylkin  [2]  on  using  wavelets  to  compress  or  sparsify 
operators.  These  results  motivate  our  use  of  wavelet  bases  to  do  multiscale  control  in  that  we 
exploit  wavelet  compression  properties  to  sparsify  our  control  transfer  functions.  This  allows 
us  to  implement  our  control  law  in  multiscale  bases  in  which  basis  coefficients  are  modified 
in  a  decoupled  fashion. 


One  of  the  main  ideas  in  using  wavelets  to  compress  operators  is  to  take  advantage  of  the  fact 
that  the  wavelet  bases  designed  in  Daubechies  [7]  can  be  used  to  efficiently  represent  smooth 
operators.  The  property  of  these  wavelets  that  allows  such  compression  is  one  that  is  referred 
to  as  the  property  of  “vanishing  moments.”  This  property  is  described  as  follows. 


x'^dx  =  0  for  m  =  0, 1,  ...M  —  1 


(18) 


It  is  also  true  that  the  number  of  vanishing  moments,  M,  is  proportional  to  the  support  length 
of  the  wavelet  or  the  length  of  the  discrete  filters  h{n)  and  g{n). 


The  reason  this  property  is  essential  to  compressing  smooth  operators  is  the  fact  that  when 
integrated  against  a  wavelet  basis  function,  ^(2”* a:  —  n),  satisfying  eq.(18),  a  function  f{x) 
is  equal  to  /(2“'”n)  plus  derivatives  of  /  of  order  greater  than  or  equal  to  M.  For  smooth 
functions,  these  higher-order  derivatives  are  small,  and  thus  the  wavelet  coefficients  are  small. 
What  this  means  for  operators  is  that  we  can  relate  compression  efficiency  of  the  wavelet 
basis  in  terms  of  smoothness  properties  of  the  operator.  Beylkin  [2]  considers  the  class  of 
Calderon-Zygmund  integral  operators,  which  have  the  property  of  being  smooth  away  from 
the  diagonal.  In  particular,  2D  operators  K{x,  y)  in  this  class  satisfy  the  following: 


1*^(2!, !/)l  < 

1 

(19) 

\x-y\ 

\d;^K(x,y)\  +  \d^K{x,y)\  < 

Cm 

(20) 

X  — 

From  the  fact  that  the  Mth  partials  of  K  decay  strongly  away  from  the  diagonals  along  with 
the  vanishing  moment  property  in  eq.(18),  one  can  show  that  the  wavelet  coefficients  decay  in 
a  similar  fashion  as  a  function  of  the  distance  between  coefficients,  i.e.,  the  distance  between 
the  centers  of  the  basis  functions  corresponding  to  two  coefficients. 


Our  control  transfer  functions  P  and  Q  will  of  course  not  be  as  simple  as  the  operators  in 
the  Calderon-Zygmund  class.  In  fact  they  are  frequency-dependent  and  depend  heavily  on 
boundary  conditions  and  material  properties  of  the  structure  on  which  the  actuators  and  sen¬ 
sors  reside.  Our  multiscale  approach  is  based  on  detailed  characterization  of  P  and  Q  in  terms 
of  their  smoothness  properties  and  any  possibly  existing  singularities.  We  will  tailor  our  mul¬ 
tiscale  bases  to  take  advantage  of  the  smoothness  of  these  transfer  functions.  The  key  issue 
is  in  trading  off  complexity  and  regularity  of  the  transform,  for  example  as  illustrated  in  the 
lattice  of  Figure  5,  and  the  efficiency  realized  by  such  a  transform,  with  the  complexity  of  the 


19 


transfer  functions  as  measured  by  their  degree  of  smoothness  as  a  function  of  position  (for 
example  position  of  actuators  and  sensors  on  the  structure). 


4.4  Designing  Neutralization  Path  P 

The  theoretical  basis  for  making  our  neutralization  transfer  function  P  sparse  in  the  wavelet 
domain  is  based  on  viewing  this  transfer  function  at  each  frequency  as  a  spatial.  Green’s 
function  operator  that  maps  inputs  (i.e.,  actuator  inputs  in  some  spatial  domain),  into  outputs 
(i.e.,  sensor  outputs  in  that  domain).  A  simple  example  of  such  an  operator  is  the  mapping 
from  normal  force  inputs  along  a  beam  to  normal  velocities  measured  along  the  beam.  Our 
analysis  applies  to  general  operators  along  surfaces;  thus,  in  particular,  we  can  treat  the  case 
of  3D  structures  with  transducers  embedded  on  their  surfaces. 

In  analyzing  the  structure  of  our  transformed  Green’s  function  operator,  K,  we  have  deter¬ 
mined  that  the  following  key  properties  of  K  can  be  used  to  compress  K  using  wavelet  trans¬ 
forms.  Note  that  K{x,y)  is  notation  for  a  general  input/output  Greens’s  function  operator 
where  y  is  variable  representing  actuator  location,  while  x  is  variable  representing  sensor  lo¬ 
cation. 

1.  K{x,  y)  is  smooth  forx^y  (i.e.,  all  its  derivatives  are  continuous). 

2.  K{x,  y)  has  a  discontinuous  Lth  partial  derivative  in  both  x  and  y  evaluated  at  the  point 
x  =  y. 

From  these  properties  we  can  determine  the  structure  of  the  transformed  K.  The  transformed 
K  maps  wavelet  transforms  of  the  inputs  to  wavelet  transforms  of  the  outputs.  At  any  partic¬ 
ular  scale  a  wavelet  coefficient  corresponds  to  the  inner  product  of  a  wavelet  function  with  a 
set  of  inputs  (outputs).  At  that  scale,  the  wavelet  coefficient  corresponding  to  a  set  of  outputs 
is  a  function  of  the  wavelet  coefficient  of  the  inputs  at  that  scale  at  locations  coincident  with 
the  outputs  and  the  wavelet  coefficients  corresponding  to  outputs  at  neighboring  locations. 
The  output  wavelet  coefficient  is  also  a  function  of  input  wavelet  coefficients  at  other  scales, 
but  again  the  dependence  involves  only  input  locations  neighboring  the  outputs.  This  local 
dependence  in  scale  and  position  characterizes  the  sparseness  of  Green’s  function  operators  in 
the  wavelet  domain,  and  is  the  basis  for  our  development  of  sparsely  interconnected  neutral¬ 
ization  filters.  A  more  precise  statement  and  more  detailed  analysis  of  the  dependence  of  the 
transformed  operator  on  scale  and  position  is  found  in  [3]. 

Figures  7  and  8  are  plots  of  the  magnitude  of  the  Green’s  function  at  1000  Hz  and  5000  Hz, 
respectively,  where  x  and  y  are  uniformly  sampled  with  128  points  over  the  interval  [0, 1]. 

Figures  9  and  10  are  plots  of  the  magnitude  of  the  wavelet  transform  of  the  Green’s  function  at 
1000  Hz  and  5000  Hz,  respectively.  Thus,  these  figures  depict  the  operator  mapping  wavelet 
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Gr66ns  function:  1000  Hz 


Figure  9:  Intensity  image  of  magnitude  of  wavelet-transformed  Green’s  function  at  1000  Hz 
(intensity  plotted  in  dB) 

coefficients  of  the  forcing  inputs  into  wavelet  coefficients  of  the  sensing  outputs.  Note  that  the 
intensities  in  these  plots  are  plotted  in  a  dB  scale  due  to  the  fact  that  coarse-scale  coefficients 
dominate  in  energy  over  many  of  the  finer-scale  coefficients. 

In  Figure  1 1  we  show  the  coefficients  of  the  wavelet-transformed  Green’s  function  at  1000  Hz 
greater  than  le  -  5  times  the  norm  of  the  Green’s  function.  Note  that  the  set  of  remaining  co¬ 
efficients  is  small,  making  the  wavelet-transformed  operator  extremely  sparse.  Similar  results 
are  obtained  for  the  Green’s  function  at  5000  Hz,  as  shown  in  Figure  12. 


4.5  Designing  Feedforward  Path  Q 

In  this  section  we  point  out  a  significant  difference  between  designing  a  multiscale  version 
of  P  and  designing  a  multiscale  version  of  Q.  In  the  case  of  the  transfer  function  P,  we 
are  simply  finding  the  best  approximation  of  the  true  plant  transfer  function  characteristics 
in  a  multiscale  basis  via  the  transformations  described  in  the  previous  section.  In  the  case 
of  the  transfer  function  Q,  however,  we  actually  design  Q  to  optimize  our  control  objective. 
Thus,  we  have  the  freedom  to  implement  Q  with  arbitrary  sparsity  and  structure,  for  example, 
trading  off  performance  for  increased  sparsity.  We  have  investigated  two  different  sparse, 
regular  structures  for  designing  and  implementing  Q.  The  first  is  the  design  of  Q  in  the  same 
multiscale  basis  as  used  to  approximate  P,  while  the  second  is  to  design  a  “banded”  Q  in  the 
spatial  domain  that  leads  to  implementing  a  locally  interconnected  feedforward  path. 
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Figure  10:  Intensity  image  of  magnitude  of  wavelet-transformed  Green’s  function  at  5000  Hz 
(intensity  plotted  in  dB) 


Figure  11:  Intensity  image  of  magnitude  of  wavelet-transformed  Green’s  function  at  1000  Hz: 
coefficients  greater  than  le-5  times  norm  of  K 


Wavelet  ttansformed  Greens  function;  magnitudes  >  1e-5  IGI 


Figure  12:  Intensity  image  of  magnitude  of  wavelet-transformed  Green’s  function  at  5000  Hz: 
coefficients  greater  than  le-5  times  norm  of  K 

4.5.1  Optimizing  Q  in  Multiscale  Bases 

We  can  rewrite  the  optimization  problem  in  eqs.(5,7)  as 

"S  \\T  +  AQS\\  (21) 

WDQWco  <  1 

where 

A  =  AYi  (22) 

S  =  Yr^S  (23) 

and  D  is  not  necessarily  related  to  D  via  a  transformation  with  Yi.  Note  that  eq.(21)  and 
eqs.(5,7)  are  equivalent  due  to  the  orthonormality  of  the  wavelet  transform  represented  by  Yi. 
In  order  to  redesign  Q  =  Y^^QYr  to  be  sparse  in  its  multiscale  representation,  we  constrain 
Q  to  have  local  interconnections  by  zeroing  out  the  appropriate  entries  of  its  matrix  represen¬ 
tation.  We  then  solve  eq.(21)  for  a  sparse  version  of  Q,  whose  interconnections  are  local. 


4.5.2  Banded  (Locally  Interconnectd)  Q 


We  define  a  banded  Q  to  be  one  whose  matrix  representation  at  each  frequency  is  banded, 
i.e..  the  matrix  is  non-zero  at  entries  some  fixed  number  of  entries  to  the  right  and  left  of 


sensor  outputs 


disturbance 


Figure  13:  Disturbance  rejection  block  diagram 


the  diagonal  and  zero  everywhere  else.  We  show  in  this  section  that  the  optimal  Q  for  a 
general  disturbance  rejection  problem  is  indeed  a  banded  operator  in  the  spatial  domain  at 
each  frequency.  We  also  show  that,  using  our  wavelet  parameterization  for  the  plant,  the 
optimal  Q-parameter  is  in  fact  “finger-like”  in  precisely  the  same  way  as  the  plant  is  in  the 
wavelet  domain.  The  thrust  of  the  development  is  based  on  finite  difference  approximations 
to  the  partial  differential  equation  (PDF)  associated  with  the  physics  of  the  plant. 

The  disturbance  rejection  problem  we  are  considering  is  the  following,  where  we  assume  P  is 
a  stable  matrix  transfer  function: 


m\n\\P{I-KP)-%  (24) 

where  K  e  H2.  The  block  diagram  for  this  problem  is  shown  in  Figure  13,  where  P  is  the 
transfer  function  from  actuator  inputs  to  sensor  outputs. 


In  Q-parameterized  form,  where  the  we  assume  the  neutralization  transfer  function  is  equal  to 
P,  the  problem  becomes 


mm||F  +  P<3P||, 


(25) 


for  Q  €  H2. 


We  show  in  the  following  that  for  systems  governed  by  linear  PDFs,  the  optimal  transfer 
function,  (J,  is  a  banded  matrix  transfer  function.  Furthermore,  the  band  is  purely  a  function 
of  the  differential  operator  associated  with  the  PDF.  We  then  show  that  banded  matrix  transfer 
functions  in  the  wavelet  domain  have  the  same  multiband  or  “finger-like”  structure  as  P  does. 
Thus,  in  the  wavelet  domain  Q  has  the  same  basic  structure  as  P. 


4.5.3  Continuous  Solutions 

Since  P  represents  the  coupling  from  actuator  inputs  to  sensor  outputs,  where  the  sensors  and 
actuators  are  located,  e.g.,  on  the  surface  of  a  structure,  we  model  P  as  a  sampled  version  of 
a  frequency-dependent  Green’s  function  integral  operator  that  maps  physical  input  variables 
on  the  surface  (e.g.,  normal  forces)  to  physical  output  variables  on  the  surface  (e.g.,  normal 
velocities).  This  Green’s  function  operator  is  associated  with  the  partial  differential  operator 
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that  contains  partial  derivatives  with  respect  to  the  space  and  time  variables  x  and  t, 
respectively.  For  example,  the  following  partial  differential  operator  corresponds  to  the  case 
of  the  Euler-Bemoulli  beam. 


$ 


x,t 


d  d 


(26) 


where  the  Cj’s  are  functions  of  x  and  t.  For  the  analysis  in  this  report,  we  consider  the  class  of 
linear  PDFs.  To  make  the  specification  of  the  PDF  well-posed  we  must  also  include  boundary 
conditions.  To  do  this  we  specify  the  initial  conditions  of  the  system  as  well  as  the  surface 
boundary  conditions. 


We  define  G{x,t,x',T)  to  be  the  Green's  function  kernel}  The  Green’s  function  integral 
operator  maps  an  input  function  /(a:,  t)  (e.g.,  the  applied  force  on  a  structure)  into  an  output 
function  «(a;,  t)  (e.g.,  the  resulting  normal  velocity  on  a  structure  due  to  the  forcing  function) 
as  follows: 

I  IG{x,t,x',T)f{x\T)dx'dT  =  u{x,t)  (27) 

The  partial  differential  operator  is  the  inverse  of  the  Green’s  function  integral  operator  in  the 
following  sense,  which  defines  the  Green’s  function. 


^x',tG{x,  t,  x',  r)  =  5{x  —  x')5{t  —  t)  (28) 

Note  that  G{x,  t,  x\  r)  can  be  viewed  as  the  causal  solution  to  an  initial  value  problem  with 
forcing  that  is  a  spatial  impulse  at  x'. 


Recall  that  Q  is  the  solution  to  the  optimization  problem  in  eq.(24).  In  the  continuous  domain, 
eq.(28)  shows  that  the  operator  $a:,<  is  in  fact  equal  to  the  optimal  solution  Q,  since  it  sets  the 
cost  function  equal  to  zero  and  is  a  causal  operator.  We  now  show  that  the  discretization  of 
the  continuous  differential  operator  is  in  fact  a  banded  matrix  with  respect  to  the  spatial 
variable  x.  Thus,  at  each  point  in  time  (or  frequency  if  represented  in  the  frequency  domain), 
the  discretized  operator  $  is  a  locally  interconnected  transfer  function. 


4.5.4  Discrete  Solutions 

For  practical  systems  we  must  model  P  as  a  sampled  version  of  a  frequency-dependent  Green’s 
function  integral  operator.  Since  the  optimal  Q  is  defined  in  terms  of  this  P,  we  analyze  the 
relationship  between  the  continuous  space/time  characterization  of  the  system  given  in  the 
previous  section  and  a  discrete  approximation  to  this.  We  do  this  by  using  n  finite  difference 
approximation  of  the  differential  operator,  This  approximation  preserves  the  banded 
structure  of  the  operator  and  in  particular  converts  the  differential  operator  into  a  banded 

^Whereas  in  the  mathematical  physics  literature  the  Green’s  function  is  associated  with  the  integral  over  the 
spatial  variables  for  a  fixed  frequency,  we  are  considering  the  function  to  vary  in  space  and  time,  with  the  integral 
taken  over  both  space  and  time  variables. 
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matrix  at  each  point  in  time.  We  will  now  show  that  in  the  discrete,  sampled  case,  the  control 
objective  function,  eq.(25),  is  small  using  the  banded  finite  approximation  of  By  small, 
we  mean  that  the  value  of  the  objective  function  in  eq.(25)  is  on  the  order  of  /i”,  where  h  is 
the  spacing  between  samples  and  n  is  the  order  of  the  differential  operator.  The  analysis  is 
performed  for  the  case  of  the  Euler-Bemoulli  beam  but  can  be  generalized  to  the  larger  class 
of  linear  PDEs. 

Consider  the  finite  difference  approximation  of  the  Euler-Bemoulli  beam  equation. 

du,  u(xq)  —  u{xo  —  A) 

^  (29) 

The  discrete  difference  operator  is  defined  as  follows. 

Dij  =  [00...d...00]  (30) 

d  =  i  [1  -8  28  -56  70  -56  28  -8  1]  (31) 

This  operator  is  in  fact  a  banded  matrix,  which  implies  a  locally  interconnected  transfer  func¬ 
tion  at  any  particular  frequency.  We  will  denote  this  matrix  as  11. 

If  the  discrete,  finite  difference  operator  11  is  well-posed,  it  has  an  inverse,  represented  by  the 
operator  P'.  Furthermore,  if  11  is  consistent  and  convergent  (this  is  true  for  our  finite  difference 
approximation  [14]),  then  its  inverse,  P',  applied  to  an  arbitrary  /  satisfies  the  following: 


u  =  Pf  (32) 

u'  =  P'f  (33) 

||«-u'||  =  0{h^)  (34) 

where  n  is  the  order  of  the  finite  difference  operator. 

ll(P-nilll/ll  <  OCA")  (35) 


Let  Q  =  — n.  Then, 


I|P(/  +  <3P)I|2  =  ||/’Q(P  -  f')||2 

(37) 

<  iiPoihiKf-nib 

(38) 

(39) 
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We  can  bound  ||P(?||2  as  follows. 


\\PQh  =  \\(P  -  P')Q  +  PQh 

(40) 

<  ll(P  -  nc?l|2  +  ll/'Qlh 

(41) 

<  ll(/’-nil2||QI|2  +  ||/||2 

(42) 

<  ||OI|20(A’-)  + 1 

(43) 

l|/’(/  +  0P)||2<O(fc")  +  ||Q||,O(A^*) 

(44) 

Eq.(44)  shows  that  the  use  of  Q  which  equals  the  discrete  approximation  of  the  differential 
operator  results  in  a  cost  function  that  is  arbitrarily  small  as  a  function  of  sensor/actuator 
sampling  distance  h.  Since  Q  is  in  fact  a  banded  matrix,  this  result  motivates  the  use  of  a 
banded  Q  for  design  of  the  optimal  feedforward  transfer  function. 


4.5.5  Structure  of  Banded  Matrices  in  Wavelet  Domain 


In  this  section,  we  show  that  banded  Q  designs  can  be  implemented  as  sparsely  interconnected 
transfer  functions  in  the  wavelet  domain.  We  define  our  banded  operator  Q  at  a  particular 
frequency  as 

Q{x,x')  =  0,  \x  -  x'\  >  {B  -  l)/2  (45) 

where  the  width  of  the  band  is  5,  x  is  the  actuator  input  position,  and  x'  is  the  sensor  output 
position. 


We  now  show  that  the  wavelet  transform  of  the  banded  operator  Q  gives  the  same  sparse 
“finger-like”  structure  that  our  neutralization  transfer  function  has.  Consider  the  following 
wavelet-transformed  operator  Q. 


ocTn,n,m,n  =  /  /  __  Q{x,  x')V’(2’"x  —  n)V>(2’"x'  —  n)dxdx' 

J2~^n  J2^^n 


(46) 


where  V’  has  M  -f- 1  vanishing  moments.  From  eq.(45)  and  the  fact  that  the  wavelet  function 
V>(2™x  -  n)  has  support  width  we  deduce  that  am,n,m,n  =  0  when  |x  -  x'|  > 

{B  —  l)/2  -1-  M.  Thus,  the  coefficients  of  Q  in  the  wavelet  domain  exhibit  the  same  sparse 
structure  that  is  exhibited  by  the  neutralization  transfer  function. 


4.6  Summary  of  Multiscale/multirate  Methodology 

A  major  benefit  of  doing  our  multiscale  decomposition  of  the  controller,  aside  from  simplify¬ 
ing  the  controller  interconnections,  is  the  fact  that  for  both  the  neutralization  and  feedforward 
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paths  the  fine-scale  components  will  in  general  have  shorter  impulse  responses,  while  coarse- 
scale  components  have  lower  bandwidth.  This  allows  shorter  FIR  filters  to  be  used  to  control 
fine-scale  components,  while  coarse-scale  components  can  be  controlled  at  a  slower  sampling 
rate.  We  summarize  the  steps  involved  in  our  multiscale/multirate  control  methodology  for 
designing  and  implementing  a  multiscale/multirate  controller  as  follows. 


4.6.1  Wavelet-Based  Neutralization  Path 

•  We  determine  our  model  of  the  plant  transfer  function  P  from  calibration  measurements 
using  an  FIR  model.  This  identified  model  forms  the  basis  of  our  neutralization  transfer 
function. 

•  We  choose  a  wavelet  basis  with  a  certain  number  of  vanishing  moments.  The  higher 
the  number  of  vanishing  moments,  the  better  the  compression  of  the  plant  P.  Recall, 
however,  that  the  support  width  of  the  wavelet  function  and,  consequently,  the  number 
of  taps  of  the  FIR  filters  h  and  g  to  implement  the  scale-recursive  wavelet  transform, 
is  proportional  to  the  number  of  vanishing  moments.  Thus,  there  is  an  inherent  trade¬ 
off  between  achieving  good  compression  and  maintaining  a  computationally  efficient 
transform. 

•  We  choose  a  sparse  interconnection  “mask”  in  the  wavelet  domain  that  represents  the 
interconnections  we  retain  for  the  neutralization  transfer  function.  This  mask  can  be 
determined  from  the  data  and  represents  neutralization  interconnections  in  the  wavelet 
domain. 

•  We  implement  the  neutralization  transfer  function  using  a  multichannel  FIR  model. 
Each  channel  will  be  low-pass  filtered  and  down-sampled  in  accordance  with  its  band¬ 
width,  while  its  impulse  response  length  will  also  vary  with  channel  in  accordance  with 
its  degree  of  damping. 

4.6.2  Q-Parameter  Design  Feedforward  Path 

•  We  do  a  Q-parameter  design  of  the  controller,  keeping  in  mind  the  implementation  of 
our  controller  in  terms  of  the  transfer  functions  Q  and  P  as  described  in  the  block 
diagram  in  Figure  4.  We  determine  the  optimal  transfer  function  Q  by  solving  the 
optimization  problem  in  eqs.(5,7)  either  (1)  in  the  wavelet  domain,  i.e.,  eq.(21),  while 
constraining  the  wavelet-transformed  operator  Q  to  have  limited  interconnections;  or  (2) 
using  a  banded  Q  transfer  function  and  solving  eqs.(5,7)  for  the  nonzero  coefficients. 

•  We  implement  the  feedforward  transfer  function  using  a  multichannel  FIR  model.  As 
in  the  case  of  the  neutralization  filter,  each  channel  will  be  subsampled  and  have  a  filter 
length  that  is  appropriate  for  that  channel.  In  contrast  to  the  neutralization  path,  however, 
for  the  feedforward  path  we  can  design  our  weights  to  have  the  multirate,  variable  filter 
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length  structure  by  further  constraining  the  structure  of  the  Q-parameter  when  solving 
the  optimization  problem  in  eq.(21). 
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5  SENSOR/ACTUATOR  PLACEMENT  METHODS  AND 
INTEGRATION  WITH  KIKO 


5.1  Overview 

The  KIKO  control  design  methodology  described  to  this  point  assumes  a  starting  point  of  uni¬ 
formly  spaced  arrays  of  collocated  sensors  and  actuators.  Application  of  the  wavelet  transform 
treats  the  model  of  the  system  near  the  boundaries  of  the  arrays  by  periodic  repetition  of  the 
array. 

In  practice  these  aspects  must  be  modified.  Generally  the  arrangement  of  transducers  is  not 
uniquely  dictated  by  the  application,  but  must  be  determined  as  part  of  the  design  process. 
Uniform  regularity  of  placement  will  generally  be  impossible  due  to  design  restrictions,  such 
as  the  presence  of  other  equipment,  projections  from  surfaces,  or  variations  in  materials.  Con¬ 
cern  for  cost  could  also  be  a  major  motivation  for  placing  transducers  more  densely  in  some 
regions  than  in  others.  Even  if  one  starts  from  an  array  designed  to  be  uniform,  yield  from  the 
manufacturing  process  and  failure  during  operation  will  both  result  in  irregular  arrays.  Man¬ 
ufacturing  tolerances  will  also  give  arrays  that  have  imperfections  in  collocation  of  sensors 
and  actuators.  Periodic  repetition  in  computations  with  a  finite  array  may  lead  to  artifacts  in 
the  size  of  wavelet  coefficients  (due  to  apparent  discontinuities  from  the  periodization)  and  a 
reduction  in  attainable  sparsity  of  the  compressed  representation.  In  the  same  vein,  similar 
artifices  would  be  necessary  to  treat  boundaries  that  are  not  rectangular,  with  similar  penal¬ 
ties.  Finally,  there  may  be  additional  benefits  in  compression  from  using  multidimensional 
wavelets  that  are  not  just  cross  products  of  one  dimensional  wavelets. 

We  have  addressed  these  issues  with  two  general  areas  of  investigation:  optimal  transducer 
selection,  and  wavelets  for  bounded  and  nonuniform  arrays. 

The  optimal  transducer  selection  work  has  resulted  in  a  computational  algorithm  to  treat  an 
appropriate  optimization  problem.  This  approach  assumes  a  predetermined  grid  of  candidate 
transducer  locations.  The  algorithm  selects  a  predetermined  number  among  these  for  use 
in  control  design.  Thus,  the  result  of  transducer  selection  will  be  an  irregularly  spaced  set 
even  if  the  set  of  candidates  is  uniform  spatially.  Our  algorithm  is  heuristic  in  order  to  avoid 
the  computational  complexity  of  an  fully  optimal  solution.  However,  our  results  do  provide 
bounds  on  suboptimality  based  on  analysis  of  the  problem  data.  Computational  examples 
illustrate  good  performance  of  the  algorithm. 

In  order  to  apply  KIKO  techniques  to  a  transducer  suite  selected  by  optimization  techniques, 
it  is  necessary  to  pursue  the  issues  of  what  to  do  with  non-uniformly  located,  possibly  non- 
collocated,  transducers  in  an  irregularly  bounded  region.  Our  approach  to  applying  wavelets 
to  these  cases  followed  three  directions.  First,  we  looked  at  direct  use  of  nonuniform  and  non- 
collocated  data  by  simply  ignoring  the  imperfections  and  proceeding  with  computation.  Next, 
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Figure  14:  Standard  feedback  control  block  diagram  with  Q-parametrized  controller. 


we  analyzed  certain  transformations  of  non-uniform  data  to  synthesize  a  uniform  array,  and  at 
the  compression  resulting  from  subsequent  wavelet  transformation.  Finally,  we  investigated 
two  types  of  generalized  wavelets:  wavelets  tailored  to  the  interval  but  using  regularly  spaced 
data,  and  so-called  second  generation  wavelets,  which  can  be  tailored  to  arbitrary  spacings  and 
domain  boundaries.  We  have  computational  results  illustrating  these  various  investigations. 

5.2  Transducer  Selection 

This  section  sununarizes  work  more  completely  described  in  [1 1],  The  problem  under  con¬ 
sideration  may  be  described  with  reference  to  Figure  14.  Because  of  the  properties  of  the 
Q-parametrization,  as  discussed  in  Section  3,  for  a  given  set  of  actuators  and  sensors  the  con¬ 
trol  design  problem  assumed  here  consists  of  choosing  the  transfer  function  Q  to  minimize. 
quadratic  cost  criterion 

<^h 

min  £  II  T(u,)  -  /l(a,)<3(u,)5(a,)  g  (47) 

^  W=Wl 

subject  to  power  constraints  on  the  actuator  commands  expressed  by 

w/i 

2  II  Q(^)S(w)  Hi  <  c  (48) 

and  the  constraint  that  Q  be  causal  and  stable.  The  precise  criteria  used  are  not  key  to  the 
developments  here,  as  will  be  seen  below.  The  u;s  in  these  expressions  are  frequencies  at 
which  we  choose  to  evaluate  the  criteria,  selected  to  cover  a  desired  bandwidth  with  sufficient 
density  to  assure  smooth  behavior  over  that  band. 
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In  the  present  discussion  we  are  concerned  with  selection  of  actuators  and  sensors,  which 
amounts  to  determination  of  the  transfer  functions  A{u})  and  S{uj).  In  the  present  approach 
we  assume  that  some  Ai{uj)  and  5'i(c(;)  are  given,  but  that  we  wish  to  reduce  the  numbers  of 
transducers  by  selecting  only  some  of  those  represented  in  Ai{u))  and  5'i(a;).  We  represent 
this  selection  by  using  the  selector  matrices  Ha  and  Hj,  which  contain  all  zeros  except  for 
a  single  one  in  each  column  of  Ilaand  each  row  of  Eg.  Ideally,  we  would  like  to  find  these 
selector  matrices  to  minimize  the  same  cost  as  above,  transforming  the  problem  into 

min  II  r(u;)  -  A(u;)na(5Hn5-S'(a;)  II2  subject  to  ^  ||  (5(w)n^5'(u;)  Hj  <  c 

(49) 

Exact  solution  of  this  problem  requires  solving  ^  ^  ^  design  problems, 

where  m  is  the  number  of  candidate  sensors,  n  the  number  of  candidate  actuators,  k  the 
number  of  desired  sensors,  and  I  the  number  of  desired  actuators.  If  the  desired  number 
of  transducers  is  flexible,  the  combinatorial  complexity  would  be  increased  by  summation  of 
candidate  values  for  k  and  I . 

As  discussed  in  [1 1],  our  approach  is  to  substitute  two  heuristic  calculation  in  order  to  avoid 
the  intractability  of  the  combinatorial  solution.  These  are 

ininmjixZ)(5'^(w),5^(a;)nJ)  and  rninmjix£)(A(u;),  A(u;)no)  (50) 

Here  D{X,  Y)  =  sin0(7?.(A’),'R(y))  is  the  sine  of  the  largest  principal  angle  between  the 
range  spaces  of  X  and  Y,  and  is  a  useful  measure  of  the  distance  between  the  ranges  of  the 
two  linear  maps  X  and  Y. 

Because  A  (and  S)  may  be  near  rank  deficient,  it  is  further  desirable  to  maintain  maximum 
numerical  rank  for  the  resulting  A(a;)na  matrix  (respectively,  n,5'(w)),  and  this  can  be  done 
by  using  the  subset  selection  algorithm  based  on  the  singular  value  decomposition  as  suggested 
in  [10].  One  benefit  of  the  subset  selection  approach  is  that  it  provides  upper  bounds  on 
the  approximation  error  in  using  a  particular  transducer  suite  in  terms  of  the  largest  residual 
singular  value.  This  provides  a  means  for  selecting  a  “good”  number  of  transducers  when 
there  is  a  gap  in  singular  values  of  of  the  A  (or  S)  matrices.  See  [11]  for  details. 

For  the  case  of  40  candidate  transducer  locations,  with  20  to  be  kept,  and  300  frequencies 
to  be  used  for  evaluation.  Table  1  indicates  the  relative  complexity  of  the  present  algorithm 
(“multidimensional  QR”)  compared  to  the  optimal  (“exhaustive  search”)  and  other  heuris¬ 
tic  approaches.  Multidimensional  QR  is  guaranteed  to  have  superior  performance  the  other 
heuristic  techniques,  and  is  roughly  comparable  in  the  complexity  of  computation. 

This  algorithm  was  used  in  the  DARPA  SPICES  program  for  transducer  selection.  [11]  has 
some  examples  of  the  SPICES  computations,  and  also  more  details  about  the  computational 
complexity. 
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Method 

Complexity 

Exhaustive  Search 

6.9el9 

Sequential  Forward  Selection 

3.8e9 

Multidimensional  QR 

9.9e7 

Best  A; 

1.7e6 

Table  1:  Example  Complexity  Comparison 


5.3  Direct  Use  of  Imperfect  Data 

As  a  first  step  in  considering  alternatives  to  compensate  for  nonuniform  or  noncollocated 
transducers,  we  computed  wavelet  transforms  of  simulated  data  for  an  Euler-Bemoulli  beam 
model  containing  the  irregularities.  Typical  results  are  shown  in  Figures  16  and  17.  For  these 
plots,  wavelet  transforms  were  computed  as  if  the  data  were  all  regularly  spaced,  when  in  fact 
the  data  contains  responses  due  to  perturbed  transducers. 


Figure  15:  8-tap  wavelet  transform  for  beam  at  1  kHz  with  regularly  spaced,  collocated  trans¬ 
ducers. 


In  Figure  16  we  see  the  results  of  perturbing  the  location  of  a  single  transducer  pair  by  one 
half  the  regular  transducer  spacing.  Note  that  one  new  vertical  and  one  new  horizontal  line 
appear  at  each  scale  in  the  transformed  data.  This  indicates  that  the  transfer  functions  from  the 
displaced  actuator  to  all  the  sensors  have  acquired  an  apparent  discontinuity  in  a  derivative,  as 
have  the  transfer  functions  from  all  the  actuators  to  the  displaced  sensor.  This  is  as  one  would 
expect:  smoothness  at  a  point  would  generally  be  disrupted  by  an  arbitrary  jump  in  a  function 
value  at  that  point. 
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Figure  16:  8-tap  wavelet  transform  for  beam  at  1  kHz  with  two  transducer  pairs  displaced 
from  evenly  spaced  position  by  0.5  the  regular  spacing. 

In  Figure  17  we  see  a  similar  result  for  the  case  of  perturbing  a  sensor  location  but  not  perturb¬ 
ing  the  location  of  the  corresponding  actuator.  Here  we  see  the  apparent  discontinuity  in  the 
transfer  functions  from  all  actuators  to  the  perturbed  sensor,  and  of  course  all  other  transfer 
functions  remain  unaffected. 

If  the  true  positions  of  transducers  are  known,  then  we  can  attempt  to  compensate  for  the  as¬ 
sumption  of  uniformly  spaced  data  in  the  wavelet  transform.  We  examine  two  approaches.  In 
the  next  section,  we  look  at  interpolation  of  irregularly  spaced  data  to  the  uniformly  spaced 
case.  In  the  following  section  we  discuss  the  use  of  wavelets  tailored  to  nonuniform  data.  In 
both  cases  the  true  positions  of  the  transducers  must  be  known,  either  by  design  or  experimen¬ 
tal  calibration. 


5.4  Transformation  of  Nonuniform  Data  for  Wavelet  Calculations 

This  section  describes  our  initial  approach  to  the  problem  of  applying  KIKO  control  technol¬ 
ogy  to  a  set  of  nonuniformly  spaced  transducers.  The  KIKO  method  of  representing  MIMO 
transfer  functions  via  lossy  compression  has  been  shown  to  yield  adequately  faithful  plant 
models  with  limited  interconnections.  By  “limited  interconnections”  we  mean  that  not  every 
output  depends  on  every  input  at  all  scales.  Further,  because  of  the  apparent  properties  of  the 
wavelet  transform,  for  certain  classes  of  plants  the  limitation  of  interconnections  is  regular  and 
predictable,  suggesting  a  conceptually  simple,  repeatable  control  architecture. 
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Figure  17:  8-tap  wavelet  transform  for  beam  at  1  kHz  with  two  sensors  displaced  from  posi¬ 
tions  collocated  with  actuators  by  0.5  the  regular  spacing. 

5.4.1  Problem  Statement 

The  wavelet-based  compression  strategy  that  KIKO  employs  takes  advantage  of  geometric 
properties  of  the  Green’s  function  describing  the  plant  to  be  modeled.  It  uses  particularly 
apt,  fast,  frequency-independent  rotations,  so  that  the  transformed  matrix  transfer  function  is 
“sparse.”  For  regularly  spaced  sensor-actuator  pairs  this  geometric  dependence  produces  the 
familiar  finger  plots  that  show  important  connections  between  input  and  output  at  a  particular 
frequency  (Section  4).  In  this  part  we  will  build  upon  these  results,  borrowing  heavily  from 
the  notation  and  hooking  into  the  machinery  for  computing  the  wavelet  coefficient  bound. 

Now  consider  the  case  where  a  sensor-actuator  pair  is  displaced  from  its  regularly  spaced 
position.  One  may  now  ask  what  happens  to  compression  and  the  finger  plot?  It  is  easy  to  see 
that  the  compression  level  is  different,  and  that  the  finger  plot  changes  qualitatively.  What  is 
less  clear  is  whether  compression  increases  or  decreases,  and  by  how  much.  The  goal  of  this 
note  is  to  explore  how  we  can  quantify  this  difference. 

The  above  paragraph  outlines  the  analysis  problem;  we  have  not  yet  attacked  the  synthesis 
counterpart.  The  main  synthesis  problem  to  be  addressed  is  the  determination  of  sensor- 
actuator  locations  so  that  compression  is  maximized,  where  the  metric  of  compression  is 
ill-defined  but  presumably  includes  some  frequency  band.  At  a  single  frequency,  we  have 
some  preliminary  results  based  on  an  extension  of  the  technique  presented  here. 

A  second  caveat  is  less  severe  but  still  worth  noting.  Analysis  of  nonuniform  compression, 
as  attempted  here,  is  different  from  reconstruction  of  a  signal  from  nonuniform  samples. 
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Wavelet-based  compression  looks  to  decompose  a  matrix  of  transfer  functions,  obtained  via 
discrete  sampling  of  an  underlying  continuous  functional.  At  no  point  do  we  reconstruct  the 
functional  from  samples. 

Although  there  is  a  broad  body  of  literature  concerning  the  reconstruction  of  functions  in  var¬ 
ious  spaces  from  their  regular  or  irregular  samples,  that  work  is  not  necessarily  applicable  to 
our  control  problem.  As  we  have  seen  in  various  examples,  the  ability  to  do  vibration  control, 
for  instance,  depends  not  on  the  faithful  reconstruction  of  arbitrary  functions  (as  limited  to  be 
in  some  space)  from  their  samples,  but  merely  with  the  correlation  of  samples  with  perfor¬ 
mance  measures.  Thus,  we  obtain  the  rule-of-thumb  in  acoustics  problems  that  the  number 
of  actuators  must  equal  or  exceed  the  rank  of  the  disturbance-to-performance  transfer  matrix, 
regardless  of  the  spatial  complexity  of  the  underlying  functional. 

All  the  results  discussed  here  refer  to  single-dimension  problems,  such  as  a  beam  rather  than  a 
plate.  Extensions  to  higher  dimensions  are  conceptually  easy  but  hindered  by  the  availability 
of  spline  or  other  interpolating  functions  in  two  or  more  dimensions. 

The  final  caveat  we  will  mention  concerns  the  collocation  of  the  sensors  and  actuator.  The 
notion  of  sensor-actuator  pairs  has  been  mentioned,  and  implies  that  there  is  exactly  one  sensor 
for  each  actuator.  Here  we  state  that  as  an  assumption,  along  with  the  additional  assumption 
that  sensors  and  actuators  are  collocated.  The  second  assumption  is  employed  for  notational 
ease  and  is  readily  removed.  In  fact,  the  irregularity  of  sensors  may  be  different  than  that  of 
actuators.  The  first  assumption  cannot  be  relaxed  without  causing  some  unnaturalness  in  the 
wavelet  transform,  since  the  input  and  output  rotations  would  be  dissimilar.  Therefore,  this 
assumption  cannot  be  relaxed  without  a  generalization  of  the  KIKO  methodology. 


5.4.2  Outline  of  Approach 

The  work  presented  here  develops  simple  machinery  to  enable  the  application  of  results,  ob¬ 
tained  earlier  for  uniform  spacing,  to  the  nonuniform-spacing  case.  Before  presenting  the 
development,  we  give  a  brief  outline  of  the  simple  underlying  idea. 

Wavelet  analysis  is  properly  defined  for  functions,  /,  defined  on  some  compact,  connected 
subset  of  the  real  line,  5  C  M.  The  wavelet  coefficients  are  obtained  by  computing  inner 
products  of  wavelet  basis  functions  with  /.  Wavelet  transforms  are  defined  for  discretely 
sampled  functions,  fk,  by  equating  the  samples  with  the  output  of  the  continuous  wavelet 
transform  of  /  at  the  finest  scale.  Subsequent  scale  coefficients  are  computed  using  the  discrete 
wavelet  transform. 

Thus  there  is  an  assumption  of  sampling  uniformity  in  the  discrete  wavelet  transform,  since  the 
wavelets  used  to  produce  fk  are  uniformly-spaced  translates  of  one  another.  This  assumption 
of  uniformity  leads  directly  to  the  KIKO  result  of  sparse  wavelet-based  representation  for 
certain  classes  of  Green’s  functions. 
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The  finite-interval  discrete  wavelet  transform  can  be  viewed  as  a  linear  transformation  from 
to  M*-.  If  the  sampling  is  nonuniform,  this  linear  transformation  is  still  well  defined, 
even  if  sparse  representations  should  no  longer  be  expected.  What  level  of  compression  will 
be  obtained?  To  address  that  question,  we  return  to  the  definition  of  the  discrete  wavelet 
transform. 

The  main  idea  is  this:  although  the  samples,  fk,  are  nonuniform,  and  thus  are  not  properly 
interpreted  as  generated  by  the  continuous  wavelet  transform  of  some  /,  we  will  define  a  new 
function,  /,  for  which  the  uniform  samples  obtained  by  computing  the  continuous  wavelet 
transform  at  the  finest  scale  are  exactly  the  fk.  Additionally,  this  function  will  be  constrained 
to  belong  to  the  same  important  spaces  as  /  (this  is  basically  a  matter  of  ensuring  sufficient 
differentiability).  We  are  now  done:  the  new  function,  /,  will  be  analyzed  using  existing 
machinery  to  determine  the  compression  level  obtained. 

That  is,  we  see  G  as  the  Green’s  function  for  an  imaginary  material  whose  regular  samples 
by  design  equal  the  irregular  samples  of  G.  Because  the  relationship  between  G  and  G  is 
known,  the  functional  form  of  G  is  also  available,  and  we  use  that  to  compute  bounds  on  the 
compression  of  G.  Finally,  since  the  wavelet  transforms  for  G  and  for  G  are  identical,  we 
obtain  bounds  for  the  compression  for  G  with  irregular  samples  Gk- 


5.4.3  Notation 

Let  A  be  a  compact,  connected  subset  of  the  real  line  which  supports  a  Green’s  function  of 
interest.  Let  r,s  e  X,  and  denote  the  Green’s  function  as  G(r,  s).  The  letter  r  stands  for 
receiver,  and  s  stands  for  source.  We  assume  we  know  a  functional  representation  for  G,  so 
that  we  may  compute  derivatives. 

Let  the  location  of  sensor-actuator  pairs  be  given  by  the  points  Xk.  They  are  not  uniformly 
spaced,  but  are,  without  loss  of  generality,  ordered  so  that  Xk-i  <  Xk.  Let  us  also  assume  that 
there  are  a  finite  number  of  transducer  locations,  1  <k  <  K,  and  the  region  of  support  for  G 
is  a  finite  interval  X  =  [A,  B].  For  simplicity  we  include  the  endpoints,  though  this  may  be 
inappropriate  for  certain  boundary  conditions.  The  extension  to  (semi)  infinite  and/or  (half) 
open  regions  of  support  is  straightforward. 

Define  a  second  compact,  connected  set,  X,  which  for  the  finite  case  above  is  the  interval 
This  is  the  set  on  which  we  will  construct  the  imaginary  Green’s  function,  G.  Thus 
we  define  r,s  G  A,  and  denote  the  Green’s  function  as  (?(r,  s).  The  transducer  locations  are 
just  the  first  K  integers:  Xk  =  k. 
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Figure  18;  Dlustration  of  example  warping  function  in  one  dimension,  g,  mapping  the  irregu¬ 
larly  sampled  space,  X,  to  the  regularly  sampled  space,  X. 

5.4.4  The  Warping  Function 

To  define  the  transform  of  G  to  a  space  of  evenly  spaced  sampling  points,  we  imagine  a  one-to- 
one  and  onto  warping  function,  g  :  X  X.  (This  development  grew  out  of  reconstruction 
theory  for  nonuniform  sampling  theory  [4, 5],  and  its  application  to  nonuniform  sampling  for 
wavelets  is  believed  novel.)  Since  X  and  X  are  intervals,  we  have  immediately  that  g{A)  =  | 
and  g(B)  =  K  +  ^.  We  also  impose  the  constraints  g{xk)  =  k,  which  implements  our 
requirement  that  the  transducer  locations  in  the  imaginary  set  are  evenly  spaced.  An  example 
of  such  a  function  is  shown  in  Figure  18. 


Since  the  function  g  is  invertible,  we  write  'y  =  g~^:X-^X.So  far,  7  is  not  well  defined. 
That  is,  there  are  an  infinity  of  possible  7’s  which  meet  the  criteria  we  have  imposed.  Let 
us  agree  to  retain  the  remaining  degrees  of  freedom  for  the  time  being,  and  investigate  the 
properties  of  7  that  can  be  deduced  at  this  point. 


5.4.5  Building  a  Useful  Imaginary  Green’s  Function 

Let  us  define  the  warped  Green’s  function  on  X  as 

(51) 

A 

Note  that,  once  7  is  specified,  G  is  well  defined.  Thus  an  imaginary  material  with  a  Green’s 
function  has  been  constructed,  with  imaginary  uniformly  spaced  samples.  It  is  now  appro¬ 
priate  to  employ  the  analysis  of  the  Green’s  function  and  its  spatial  derivatives  to  bound  the 
level  of  compression,  because  the  assumption  of  sample  uniformity  has  been  met,  enabling 
a  discretely  sampled  function,  fk,  to  be  equated  with  the  output  of  the  continuous  wavelet 
transform  of  some  /  at  the  finest  scale. 
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In  [3],  a  bound  for  the  wavelet  coefficient,  otm,n,fh,n,  is  sought.  (Here  m  is  the  scale  index  and 
n  is  the  translation  index  for  the  “receive”  wavelet,  and  m  is  the  scale  index  while  h  is  the 
translation  index  for  the  “source”  wavelet.  See  that  paper  for  details.)  To  compute  this  bound, 
the  anchor  points  ro  =  2~’"n  and  io  =  were  used,  representing  the  edges  of  the  receive 
and  source  wavelet  supports,  Q,f  and  (Is,  respectively.  In  general,  they  may  represent  wavelets 
at  different  scales,  and  may  or  may  not  intersect.  Additionally,  r  and  s  are  test  points  inside 
(If  and  (1$;  these  points  will  be  varied  as  the  Green’s  function  is  integrated  against  the  wavelet 
basis  functions  to  form  the  bound.  In  summary. 


(If 

= 

[2-’”n,2-’”(n-M)] 

= 

[ro,ro  +  2-’"] 

(52) 

(Is 

=z 

[2-’^n,2-’^(n-l-l)] 

= 

[Jo;  So  +  2  ’"j 

(53) 

f 

G 

(If 

S 

G 

(Is 

The  essential  component  to  the  bound  on  compression  is  a  term, 

s)  =  ^sup^  -  ^o),  So  +  t{s  -  So))  j  (54) 

Here  M  is  the  number  of  vanishing  moments  possessed  by  the  wavelet  basis  functions,  and  i 
varies  from  0  to  2M.  This  expression  is  correct  when  (If  D  Os  is  empty;  a  slightly  different 
expression  is  required  when  the  intersection  is  nonempty: 

Recall  that  the  class  of  Green’s  functions  presented  in  [3]  are  discontinuous  in  their  deriva¬ 
tive  at  r  =  s.  The  central  difference  in  the  two  cases  is  that  when  r  =  s,  as  can  only  happen 
in  the  intersecting  case,  the  wavelet  transform  will  uncover  the  discontinuity  in  the  T*^-order 
partial  derivative  of  G,  if  2M  >  T  —  1. 

Computation  of  5)  in  (54)  is  straightforward  once  7  is  chosen.  For  notational  conve¬ 
nience,  define  the  following  identical  functions  to  serve  as  placeholders:  p(-)  =  7(.)  when  the 

argument  is  r;  <t(-)  =  7(-)  when  the  argument  is  s.  Also  define  r(x;  xo)  =  xq  -|-f(x  —  xo).  Then 
the  form  of  the  arguments  of  G  in  (54)  becomes 

G{ro  +  t(f  -  fo),  5o  +  t{s  -  Jo))  =  G{p{ro  +  t{r  -  rO)),  cr(.So  +  t{s  -  So))) 

=  G(p(r(f;fo)),a(T(s;so))) 

Partial  derivatives  of  this  term  may  be  computed  directly  using  the  chain  rule.  For  example, 

Q 

^G(p(r(f;fo)),a(r(J;so)))  =  ta'G^{p,(T) 
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^G(/9(r(r;  ro)),  <T(r(i;  5o)))  =  t^(T"Ge{p,  a)  +  a) 

-^G{p{T{f-,  ro)),  o-(r(s;  sq)))  =  t^cr"'G„{p,  a)  +  3t^a'a"G^^{p,  a)  +  t^{cr'fG^^„{p,  a) 
^^G(/!)(r(r;ro)),a(r(s;5o)))  =  t^p  a'Gp^{p,a) 

where  the  arguments  r(r;  ro)  and  r(s;  so)  have  been  suppressed.  Subsequent  derivatives  may 
be  computed  directly  by  repeated  use  of  the  chain  rule  and  the  product  rule.  The  remainder 
of  the  wavelet  coefficient  bound  computation  is  identical  to  that  given  in  the  corrected  version 
of  [3]. 


5.4.6  Bounds  Comparison 

The  comparison  of  bounds  is  unfortunately  complicated  by  the  lack  of  intuition  about  what 
factors  contribute  significantly.  For  instance,  the  term  i5,  (f ,  s)  must  be  evaluated,  either  inside 
the  integral  in  (43)  of  [3],  or  separately  to  compute  a  supremum  over  (f ,  s)  6  fir  x  which 
can  then  be  brought  outside  the  integral.  The  same  applies  for  5,  (r,  i)  in  (55)  of  [3].  We  have 
a  poor  understanding  of  the  relative  magnitudes  of  these  two  terms. 

In  particular,  there  is  an  assumption  that  Bi  is  not  so  much  larger  than  B,  that  it  contributes  to 
the  bounds  on  size  of  the  wavelet  coefficients  for  nonintersecting  and  intersecting  cases.  If  this 
assumption  were  violated,  it  could  increase  the  number  of  vanishing  moments,  M,  required  to 
get  good  compression.  The  current  mathematical  machinery  does  not  address  this  question. 

Nevertheless,  the  example  computations  above  show  that  the  computation  of  the  5,’s  or  Bi’s 
is  changed  by  nonuniform  spacing  only  through  the  introduction  of  derivatives  of  this  warping 
function.  The  warping  function  really  adds  nothing  that  essentially  complicates  computation 
of  the  bound.  Recall  that  there  were  an  infinity  of  warping  functions  that  interpolate  the 
transducer  locations  in  Figure  18.  Good  choices  for  7,  which  don’t  change  either  Bi  or  Bi 
much,  would  then  have  first  derivatives  close  to  unity  and  all  other  derivatives  small  or  zero 
throughout  the  domain  of  7.  Intuitively,  one  may  think  that  these  derivatives  never  get  larger 
than  they  have  to,  while  still  growing  the  function  to  meet  the  boimdary  conditions,  so  that 
various  products  of  them  with  G  and  its  derivatives  remain  small. 


5.4.7  Issues 

Note  that  the  expression  (54)  becomes  unwieldy  if  7  has  nonvanishing  high-order  derivatives, 
because  repeated  use  of  the  product  rule  is  required.  Each  use  of  the  product  rule  produces 
more  terms  that  do  not  group  together  in  a  compact  structure  such  as  a  binomial  expansion. 
However,  if  7  is  an  affine  function,  so  that  p  and  a  are  both  affine,  then  the  computation  of  Bi 
and  Bi  are  greatly  simplified. 
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The  choice  of  warping  function  plays  a  key  role  in  our  analysis,  but  its  effect  on  the  com¬ 
pression  bounds  is  distressingly  hidden,  and  its  effect  on  the  true  compression  is  of  course 
nonexistent.  Thus  the  act  of  warping  function  choice  becomes  one  of  preserving  the  tightness 
of  the  bounds,  if  indeed  that  is  possible. 

The  warping  function  can  be  considered  to  interpolate  among  the  {xk,  Xk  =  k)  pairs.  The  plot 
in  Figure  18  used  cubic  spline  interpolation  as  an  example.  One  could  also  consider  the  use  of 
linear  interpolation,  which  offers  the  attractive  quality  of  vastly  simplifying  the  computation  of 
the  bounds  for  almost  all  values  of  i.  However,  a  linear  interpolation  introduces  discontinuities 
in  the  first  derivative  at  the  knot  points.  It  is  not  clear  how  to  handle  these  within  the  supremum 
framework  we  have  developed;  they  clearly  cannot  be  ignored,  since  a  smooth  approximation 
to  such  a  discontinuity  would  contribute  large  factors  corresponding  to  higher  derivatives.  On 
the  other  hand,  they  are  not  real,  and  so  should  not  be  expected  to  hurt  compression  the  way 
discontinuities  in  the  Green’s  function  do. 

One  view,  requiring  further  exploration,  is  that  the  use  of  linear  interpolation  may  lead  to 
a  good  approximation  of  the  strict  bound.  This  approximation  may  prove  useful  since  it 
would  approximately  quantify  the  loss  of  compression  as  a  function  of  transducer  spacing. 
The  formula  for  computing  the  bound,  which  is  undefined  at  the  knots,  will  need  to  repaired. 


5.5  Generalized  Wavelets 

The  basic  theory  of  wavelets  applies  to  functions  in  L^(m),  that  is,  square-integrable  functions 
defined  on  the  real  line.  It  is  straightforward  to  extend  this  theory  to  the  higher  dimensional 
spaces  by  simply  taking  as  the  wavelet  basis  the  direct  product  of  wavelets  for  L^(M). 

However,  realistic  domains  for  posing  control  problems  will  be  bounded  in  M",  and  general¬ 
ization  in  a  satisfactory  way  to  other  spaces,  starting  with  X^([0, 1]),  is  not  so  obvious. 

Two  issues  arise  immediately  in  considering  the  application  of  wavelet  representations  to  real 
control  problems:  First,  that  we  need  to  compute  wavelet  transforms  for  functions  on  bounded 
domains,  and  second,  that  the  spatial  distribution  of  the  data  we  wish  to  transform  will  gen¬ 
erally  not  be  uniform  or  regularly  spaced.  The  first  issue  is  a  problem  for  standard  wavelet 
analysis  because  any  smooth  wavelet  basis  for  L^(m)  has  elements  that  straddle  the  bound¬ 
ary  points  0  and  1,  and  calculation  of  the  corresponding  coefficients  in  the  wavelet  expansion 
requires  knowledge  of  function  values  outside  the  interval  of  definition. 

Considering  procedures  used  to  handle  Fourier  transforms  on  finite  intervals,  there  are  two  ob¬ 
vious  approaches  to  try:  extending  the  function  by  padding  with  zeros,  and  extending  the  func¬ 
tion  by  periodic  repetition.  Following  [6],  we  note  several  difficulties  with  these  approaches. 
Both  the  zero  padding  and  periodic  extension  techniques  typically  generate  discontinuities  at 
the  boundary,  resulting  in  large  wavelet  coefficients  for  fine  scales  even  when  the  funciton  is 
otherwise  smooth  on  [0, 1].  Furthermore,  the  zero  padding  approach  uses  more  wavelets  than 
necessary,  since  some  wavelets  extend  outside  the  interval  in  their  support. 
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[6]  also  mentions  the  use  in  image  analysis  of  extension  beyond  the  boundaries  by  reflection, 
which  obviously  results  in  continuity  at  the  boundary  points.  However,  the  derivatives  at  the 
boundaries  will  generally  be  discontinuous  even  when  the  function  is  smooth,  again  resulting 
in  (possibly  unnecessarily)  large  wavelet  coefficients.  A  fourth  alternative  mentioned  in  [6]  is 
a  wavelet  basis  for  the  interval  developed  by  [13],  which  does  not  give  boundary  artifacts  as 
do  the  others.  However,  the  authors  of  [6]  object  to  the  inability  to  generalize  these  wavelets  to 
wavelet  packets,  and  numerical  ill-conditioning  involved  in  constructing  the  “edge”  wavelets 
with  support  at  the  boundary. 

[6]  introduces  a  new  wavelet  construction  that  does  not  suffer  from  the  above  difficulties,  and 
we  performed  some  experiments  with  these  wavelets  for  our  beam  and  plate  models.  ([6] 
cites  a  similar  construction  developed  independently  by  B.  Jawerth,  and  discussed  in  [1].) 
Other  (less  useful)  constructions  are  due  to  various  other  authors.  However,  all  these  have  in 
common  that  for  application  to  sampled  data  via  discrete  version  of  the  transforms,  uniformly 
spaced  data  are  required. 

More  recently,  W.  Sweldens  [15]developed  a  new  approach  to  the  construction  of  wavelets 
that  applies  not  only  to  functions  on  intervals,  but  to  general  domains  in  multiple  dimensions 
with  unevenly  space  data,  which  he  called  [16]  “second  generation”  wavelets.  We  also  exper¬ 
imented  with  these  wavelets  using  our  beam  and  plate  models. 

Below,  we  report  on  these  two  sets  of  experiments,  and  discuss  our  perceptions  of  their  com¬ 
parison,  advantages  and  shortcomings.  We  wish  to  emphasize  that  these  results  are  largely 
superficial  and  preliminary,  and  we  feel  that  there  is  much  benefit  to  be  derived  from  further 
investigation  in  this  area. 


5.5.1  Comparisons  of  Reconstruction  Error  and  Sparsity 

In  numerical  experiments  using  the  wavelets  constructed  in  [6],  which  we  call  “boundary 
wavelets,”  it  is  easy  to  demonstrate  significant  advantages  in  compression  over  applying  stan¬ 
dard  wavelet  analysis  through  periodic  repetition  of  a  function  on  an  interval.  However,  this 
advantage  relies  on  the  presence  of  discontinuities  at  the  (former)  boundary  points  in  the  re¬ 
peated  function.  The  frequency  domain  beam  model  used  in  our  numerical  experiments  had 
clamped  boundary  conditions,  and  therefore  the  periodized  transfer  function  is  spatially  con¬ 
tinuous  at  the  (former)  boundary  points  up  to  the  first  derivative.  This  results  in  no  advantage 
for  the  boundary  wavelets  over  standard  wavlets,  as  can  be  seen  in  Table  2.  The  data  in  Table 
2  is  based  on  setting  to  zero  all  wavelet  coefficients  below  a  certain  threshold  fraction  of  the 
norm  of  the  wavelet  transform,  listed  as  the  “truncation  level.”  The  “reconstruction  error” 
listed  is  the  fraction  of  the  norm  of  the  original  transfer  function  data. 

Table  2  also  indicates,  at  least  for  this  example,  that  the  lifted  wavelets  can  provide  a  factor 
of  2  more  sparsity  for  a  given  level  of  reconstruction  error.  It  is  difficult  to  make  direct  com¬ 
parisons,  because  of  the  number  of  variables  involved.  For  a  given  wavelet  order,  the  lifted 
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Standard  Wavelets 

Boundary  Wavelets 

Lifted  Wavelets 

Wavelet 

Order 

Tmncation 

Level 

Density 

Reconstrac- 
tion  Error 

Density 

Reconstmc- 

tion  Error 

Density 

Reconstrac- 
tion  Error 

4 

-50  dB 

24.7% 

-39.2  dB 

24.7% 

-39.0  dB 

15.2% 

-50.8  dB 

4 

-40  dB 

9.9% 

-29.7  dB 

9.6% 

-31.1  dB 

5.0% 

-38.0  dB 

6 

11.1% 

-37.9  dB 

10.8% 

-36.2  dB 

-44.0  dB 

6 

4.5% 

-30.0  dB 

4.4% 

-30.5  dB 

-30.4  dB 

Table  2:  Compression  results  based  on  128  transducer  pairs  and  4  levels  of  wavelets  computed. 


wavelet  algorithm  can  require  better  than  a  factor  of  2  less  computation.  (Since  lifted  wavelets 
are  not  orthonormal,  there  are  two  order  numbers  to  specify  -  the  number  of  vanishing  mo¬ 
ments  and  the  number  dual  vanishing  moments.  For  comparison’s  sake  here  we  have  chosen 
wavelets  with  these  two  numbers  equal.  In  practice,  however,  this  afford  an  additional  degree 
of  freedom  in  performing  the  tradeoff  between  sparsity  and  reconstruction  error.)  Different 
wavelets  may  be  chosen  for  all  three  cases,  and  performance  may  depend  on  those  choices. 
It  is  not  clear  from  this  example  whether  the  better  performance  for  the  lifted  wavelets  is  due 
to  the  treatment  of  edge  effects  or  the  difference  in  the  interior  wavelets.  Ultimately,  the  best 
basis  of  comparison  would  be  the  amount  of  computation  required  to  achieve  desired  levels 
of  sparsity  and  reconstruction  error.  We  have  not  gone  to  the  required  level  of  detail  for  that 
analysis. 


5.6  Conclusions  and  Recommendations  for  Further  Work 


The  transducer  selection  algorithm  we  have  described  above  and  in  the  appended  paper  [11] 
seeks  to  pick  transducers  that  give  the  best  numerical  range  in  sensing  and  actuation  for  fixed 
number  of  transducers.  Such  a  criterion  makes  sense  from  the  standpoint  of  minimizing  gains 
to  enhance  robustness,  and  limiting  the  sensitivity  of  performance  to  noise  or  random  errors 
in  measurement.  The  algorithm  will  result  in  nonuniform  transducer  locations,  and  of  course 
real  transducer  arrays  will  always  be  spatially  bounded.  We  discussed  two  means  of  treating 
the  nonuniformities,  interpolation  and  second  generation  wavelets,  and  two  means  of  treat¬ 
ing  boundedness,  boundary  wavelets  and  second  generation  wavlets.  It  is  clear  that  second 
generation  wavelets  merit  further  investigation. 

The  definition  of  second  generation  wavelets  will  depend  on  the  transducer  locations,  so  de¬ 
velopment  of  the  algorithms  for  performing  transforms  and  the  architecture  for  information 
flow  must  take  into  account  the  sequence  and  assumptions  behind  such  data.  For  example,  it 
is  likely  that  quite  accurate  transducer  locations  will  be  known  prior  to  fabrication,  allowing 
algorithm  structure  and  data  network  architecture  to  be  determined  prior  to  fabrication,  but 
that  it  may  be  necessary  to  perform  calibration  on  the  as-built  system  to  determine  the  variety 
of  coefficients  needed  to  complete  implementation.  Of  course,  depending  on  variability,  it 
may  be  possible  to  calibration  once  for  an  entire  production  line. 
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There  is  another  possible  use  of  transducer  selection  or  placement  which  we  have  not  inves¬ 
tigated,  that  of  selecting  transducers  to  optimize  compression.  This  makes  sense  only  when 
there  are  a  large  number  of  transducers,  and  where  much  of  the  data  is  redundant  in  an  ap¬ 
propriate  sense.  (Otherwise,  transducer  location  is  essentially  dictated  by  what  one  needs  to 
know  and  to  affect,  and  algorithm  presented  above  is  very  relevant.) 

This  reasoning  might  seem  to  lead  to  dispensing  with  large  numbers  of  transducers  and  simply 
using  the  minimally  few  devices  needed  as  determined  via  our  selection  algorithm:  this  is 
another  form  of  compression  -  done  by  installing  only  those  transducers  a  model  says  are 
needed.  However,  there  is  a  big  difference  between  compression  based  on  just  placing  devices 
where  a  model  says,  and  compression  using  the  real  smoothness  properties  of  many  physical 
measurements:  the  former  is  much  more  susceptible  to  the  effects  of  uncertainty  and  much 
less  adaptable  to  changes  in  the  system.  The  availability  of  proliferated  sensing  and  actuation 
opens  up  an  opportunity  for  substituting  real  time  measurements  for  models  that  are  dubious 
(as  in  the  case  of  complex  physical  systems). 

From  this  perspective,  it  seems  that  it  should  be  important  to  investigate  patterns  of  transducer 
arrays,  in  particular  how  to  position  devices  around  boundaries  and  physical  discontinuities, 
and  densities  of  transducer  arrays  needed  to  measure  and  affect  desired  physical  properties. 
These  are  areas  for  continuing  work. 
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6  NUMERICAL  RESULTS  AND  SIMULATIONS  FROM 
MODEL  DATA 


In  this  section  we  report  our  results  on  applying  the  KIKO  control  method  to  simulated  data 
using  two  different  numerical  models:  (1)  a  model  based  on  solutions  to  an  Euler-Bemoulli 
beam  and  (2)  a  model  based  on  modal  representation  of  a  flexible  plate. 


6.1  Euler-Bernoulli  Beam  Results 

We  used  an  Euler-Bemoulli  beam  model  to  test  our  transfer  compression  method  and  to  test 
our  overall  control  method  using  simulated  data.  For  our  example,  we  use  an  Euler-Bemoulli 
beam,  whose  equation  of  motion  is  given  by 

Wrrrr  +  =  —p{r)  T  €  [0,  1] 

where  w  is  the  dimensionless  vertical  displacement  of  the  beam,  ui  is  dimensionless  frequency 
and  p  represents  the  forcing  along  the  beam.  We  use  clamped  boundary  conditions 

ty  =  to,.  z=  0  for  r  =  0, 1. 

The  Green’s  function  in  this  case  satisfies 

G{r,  s)rrrr  +  s)  =  -5{r  -  5), 

along  with  the  boundary  conditions  above.  Solution  for  the  Green’s  function  is  possible 
in  closed  form.  An  analytical  derivation  for  the  Green’s  function  of  a  finite-length  Euler 
Bernoulli  beam  is  given  in  Appendix  A.  Figure  19  shows  the  magnitude  of  the  Green’s  func¬ 
tion  at  a  dimensionless  frequency,  u;  =  215,  where  512  equally  spaced  samples  along  the 
beam  were  used  as  collocated  forcing  and  sensing  points. 

The  result  of  applying  an  orthonormal  wavelet  transform  to  the  Green’s  function  depicted  in 
Figure  19  is  shown  in  Figure  20.  Note  that  the  coefficient  energy  is  highly  concentrated  along 
“finger-like”  patterns  in  the  matrix.  This  provides  iimnediate  visual  validation  for  our  analyti¬ 
cal  predictions  of  transfer  function  compression  by  showing  that  in  the  wavelet  representation 
the  transfer  function  coefficients  are  concentrated  into  a  sparse  pattern. 

In  Figure  21,  we  show  the  result  of  thresholding  the  wavelet  coefficients  shown  in  Figure  20 
so  that  only  coefficients  greater  than  le-8  times  the  norm  of  the  matrix  appear  in  black.  Note 
that  the  resulting  coefficient  pattern  is  indeed  extremely  sparse  and  that  these  coefficients  can 
be  used  to  approximate  the  original  Green’s  function  operator  with  a  relative  error  of  only 
le-8.  To  compare  what  our  analytical  results  in  [3]  suggested  as  important  coefficients  for 
approximation  of  the  Green’s  function  operator.  Figure  22  shows  the  predicted  finger-like 
pattern  where  overlaps  in  wavelet  function  supports  occur. 
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Figure  19:  A  plot  of  the  magnitude  of  the  Green’s  function  for  our  Euler-Bemoulli  beam  with 
clamped  ends  as  a  function  of  the  field  point,  r  and  the  source  point,  s,  and  dimensionless 
frequency,  u)  =  215.  The  scale  is  given  by  the  colorbar  on  the  right. 
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Figure  20:  Magnitude  of  wavelet  transform  of  Euler-Bemoulli  beam  Green’s  function  using 
Daubechies  8-tap  filter,  plotted  on  dB  intensity  scale 


Figure  21 


Greens  function:  1 000  Hz  x  1 0 


Figure  23:  Intensity  image  of  magnitude  of  Green’s  function  at  1000  Hz 


Further  numerical  experiments  showed  that  the  compression  property  holds  for  a  wide  range 
of  frequencies.  As  an  example,  consider  the  case  in  Figure  23,  which  shows  the  magnitude  of 
the  Green’s  function  at  a  much  lower  frequency  than  the  previous  case.  Also,  the  number  of 
samples  is  now  128  points  rather  than  512. 

Figure  24  is  a  plot  of  the  magnitude  of  the  wavelet  transform  of  the  Green’s  function.  Again, 
the  finger-like  pattern  appears.  In  Figure  25,  we  show  the  coefficients  of  the  wavelet-transformed 
Green’s  function  greater  than  le  -  5  times  the  norm  of  the  Green’s  function. 

We  now  show  results  of  controlling  the  response  of  the  Euler-Bemoulli  beam  in  the  presence 
of  external  disturbances.  These  results  show  that  wavelet  transfer  function  compression  com¬ 
bined  with  sparse  Q-parameter  design  can  be  successfully  applied  to  the  control  of  a  beam 
using  a  dense  array  of  actuators  and  sensors.  The  beam  is  taken  to  be  steel  of  length  1  m, 
thickness  2  cm  and  width  5  cm.  The  internal  damping  J7  =  .2,  while  the  Young’s  modulus  is 
equal  to  1.96ell  kg/ms^. 

We  design  an  optimal  banded  Q-parameter  using  an  FIR  parameterization  for  Q  and  solving 
eqs.(5,7)  for  the  non-zero  coefficients.  Figure  26  shows  the  performance  results  for  two  dif¬ 
ferent  disturbances:  (1)  an  impulsive  force  excitation  at  a  point  on  the  beam  and  (2)  a  random 
distribution  of  force  excitations  along  the  beam.  We  used  32  equally  sampled  points  along  the 
beam  as  collocated  sensor  and  actuator  positions.  Full  neutralization  was  used,  i.e.,  no  wavelet 
compression  was  applied  to  the  plant  transfer  function.  This  case  was  performed  to  give  us 
baseline  performance  for  comparison  purposes.  The  RMS  performance  was  19.7  dB  over  a  0 
to  10  kHz  bandwidth.  Note  that  the  bottom  two  plots  bearing  a  resemblance  to  lips  show  the 
maximum  output  response  for  the  controlled  (white)  and  uncontrolled  (red)  cases  along  the 
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Wavelet  transformed  Greens  function:  1000  Hz  (dB  intensrty  scale) 
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Figure  24:  Intensity  image  of  magnitude  of  wavelet-transformed  Green’s  function  at  1(XX)  Hz 
(intensity  plotted  in  dB) 


Wavelet  transformed  Greens  function:  magnitudes  >  1e-5  IGI 
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Figure  25:  Intensity  image  of  magnitude  of  wavelet-transformed  Green’s  function  at  KXK)  Hz: 
coefficients  greater  than  le-5  times  norm  of  K 
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length  of  the  beam,  while  the  black  bars  indicate  the  forcing  input. 

Figure  27  shows  the  performance  results  for  the  same  cases,  where  the  neutralization  transfer 
function  is  approximated  using  our  wavelet  compression  approach.  In  this  case,  the  wavelet 
sparsification  used  was  such  that  48  percent  of  the  coefficients  were  retained. 

Finally,  in  Figures  28  and  29  we  show  performance  results  in  the  presence  of  noise  for  20 
dB  and  10  dB  SNRs,  respectively,  where  the  neutralization  transfer  function  is  approximated 
again  using  the  same  48  percent  of  the  coefficients. 


6.2  Flexible  Plate  Results 

In  this  section  we  illustrate  our  control  approach  for  the  case  of  controlling  vibrations  in  a  1  m 
X  1  m  flexible  plate  with  a  16  by  16  array  of  transducers.  This  example  focuses  on  the  ability 
of  the  wavelets  to  compress  the  plant  P,  which,  in  the  case  of  our  example,  is  taken  to  be  the 
operator  that  maps  point  forces  uniformly  distributed  on  the  plate  into  normal  velocities  at  the 
same  locations.  Our  control  objective  is  to  minimize  the  vibrations  at  the  sensor  outputs  due 
to  external  disturbance  forces  by  designing  a  feedback  controller  whose  inputs  are  the  sensor 
outputs  on  the  plate  and  whose  outputs  are  the  control  signals  to  the  actuator  inputs.  This  maps 
directly  into  the  problem  of  minimizing  the  norm  of  the  transfer  function  in  eq.(l),  where  the 
transfer  functions  T,  S,  and  A  are  all  equal  to  P. 

Figure  30  is  a  plot  of  the  velocity  response  of  the  plate  due  to  a  point  force  in  the  middle 
of  the  northwest  quadrant  of  the  plate  at  250  Hz.  The  overall  transfer  function  P,  whose 
inputs  include  actuator  point  forces  sampled  uniformly  throughout  the  plate  and  whose  out¬ 
puts  include  sensor  velocity  measurements  at  points  collocated  with  the  actuators,  maps  for 
each  frequency  all  possible  actuator  inputs  to  all  possible  sensor  outputs.  Thus,  at  a  single 
frequency  we  can  think  of  P  as  a  collection  of  plots  as  in  Figure  30,  one  for  each  of  the  16 
actuator  locations.  Figure  31  depicts  this  at  250  Hz.  One  can  think  of  this  plot  as  a  16  by  16 
checkerboard,  where  each  square  shows  the  velocity  response  due  to  an  input  whose  location 
corresponds  to  the  location  of  the  square  with  respect  to  the  checkerboard.  In  Figure  31,  we 
also  show  the  wavelet-transformed  P,  which  is  obtained  using  a  4-dimensional  transform  de¬ 
rived  from  tensor  products  of  ID  transforms.  One  immediately  notices  that  the  amplitudes  of 
the  transfer  function  in  the  wavelet  domain  are  much  more  concentrated  into  a  sparse  set  of 
locations.  In  fact  if  we  threshold  the  amplitudes  to  retain  only  10  percent  of  the  channels,  we 
have  the  channels  depicted  in  Figure  32.  Using  these  channels  as  an  approximation  of  P  for 
our  neutralization  filter,  we  get  an  8  percent  error. 

We  now  show  how  the  approximation  error  in  our  wavelet-approximated  neutralization  filter 
affects  closed-loop  control  performance.  In  particular,  we  show  the  effects  of  neutralization 
error  on  a  Q-parameter  design  using  a  nearest-neighbor  structure.  For  this  structure  each  actu¬ 
ator  control  signal  input  is  a  function  of  the  sensor  outputs  at  each  of  the  eight  nearest  locations 
surrounding  the  actuator  plus  the  output  of  the  sensor  colocated  with  the  actuator.  Figure  33 
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Figure  29:  Frequency  and  time  domain  plots  of  controlled  vs.  uncontrolled  responses,  using 
banded  Q-parameter  and  wavelet  compressed  neutralization:  10  dB  SNR 
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Figure  32:  Nonzero  channels  of  wavelet-transformed  P  for  10  percent  density 
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Figure  33:  Closed-loop  performance  at  250  Hz:  (1)  perfect  neutralization  and  (2)  wavelet- 
transformed  neutralization 

shows  the  closed-loop  control  performance  at  250  Hz  using  our  designed  Q-parameter  with 
(1)  perfect  neutralization  and  (2)  wavelet-approximated  neutralization.  These  plots  are  nor¬ 
malized  to  be  compared  to  the  uncontrolled  response  shown  in  the  first  plot  of  Figure  3 1 .  The 
rms  performance  over  a  0-1550  Hz  band  for  case  (1)  is  9.9  dB,  while  for  case  (2)  it  is  9.1 
dB.  Thus,  we  only  get  a  0.8  dB  degradation  in  performance  using  a  neutralization  filter  with 
only  10  percent  of  the  original  channels.  In  other  words,  this  is  a  90  percent  reduction  in 
interconnection  complexity  with  very  little  degradation  in  performance  relative  to  using  full 
intercoimections. 

To  motivate  how  the  multirate  structure  for  the  neutralization  filter  would  work,  we  show 
the  impulse  response  of  several  channels  of  the  wavelet-approximated  neutralization  filter. 
Figure  34  shows  impulse  responses  for  (1)  several  channels  of  P  and  (2)  several  channels  of 
the  wavelet-transformed  P  corresponding  to,  from  top  to  bottom,  coarse,  medium,  and  fine 
scales,  respectively.  Note  that  the  coarser-scale  channels  exhibit  lower  bandwidth  responses, 
while  finer-scale  channels  exhibit  shorter  responses.  This  structure  forms  the  basis  of  the 
multirate  scheme  outlined  in  the  previous  section. 


6.3  Conclusions 

We  have  described  an  approach  to  controlling  large  numbers  of  sensors  and  actuators  via  a 
multiscale/multirate  processing  scheme.  By  implementing  the  controller  in  the  Q-parameterization 
we  have  identified  two  complementary  sources  of  a  processing  bottleneck.  In  the  neutraliza¬ 
tion  path  of  the  controller,  the  problem  lies  in  the  interconnections  between  actuators  and 
sensors,  while  in  the  feedforward  path  of  the  controller,  the  problem  is  in  the  interconnections 
between  sensors  and  actuators.  For  many  physical  systems,  the  neutralization  path  transfer 
function  P,  which  represents  the  coupling  between  actuators  and  sensors,  is  predominantly 
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Figure  34:  Impulse  responses:  (1)  P  and  (2)  wavelet-transformed  P 


smooth,  and  compression  can  be  obtained  by  performing  a  wavelet  transform  on  the  actuators 
and  sensors,  then  implementing  a  sparse  operator  between  the  coefficients.  Results  on  this 
were  given  for  the  both  the  case  of  a  flexible  beam  as  well  as  the  case  of  a  flexible  plate. 

The  feedforward  path  transfer  function  Q  can  also  be  designed  in  the  multiscale  domain,  giv¬ 
ing  us  the  same  sparse  interconnectivity  as  in  the  case  of  the  neutralization  transfer  function. 
For  both  P  and  Q,  the  multiscale  structure  induces  a  multirate  structure  in  the  time  domain 
that  gives  further  computational  efficiencies.  In  particular,  coarse-scale  components  can  be 
processed  at  lower  sampling  rates,  while  fine-scale  components  are  processed  at  higher  rates 
but  with  shorter  filters. 

In  summary,  the  savings  afforded  by  our  overall  multiscale/multirate  methodology  in  proces¬ 
sor  communication,  computation,  and  its  inherent  scalability,  offers  promise  for  controlling 
spatially  varying,  wideband  phenomena.  The  methodology  also  offers  a  reasonable  scheme 
for  allowing  the  controller  to  optimally  process  a  dense  array  of  transducers  for  a  given  prob¬ 
lem,  obviating  the  need  for  customizing  the  transducer  layout  to  a  particular  problem. 
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7  RESULTS  OF  EXPERIMENTAL  VALIDATION 


7.1  PSU  Experimental  Setup 


This  section  summarizes  the  data  acquisition  procedure  used  in  the  experimental  investigation 
of  the  dynamic  response  of  a  fully  clamped  panel.  The  investigation  was  conducted  in  order 
to  provide  SRI  with  experimentally  measured  structural  transfer  functions  of  a  panel  such  that 
realistic  data  is  available  for  the  design  of  a  multichannel  control  circuit.  The  measurements 
were  performed  at  the  Center  for  Acoustics  and  Vibration  Noise  Control  Laboratory  located 
at  the  Penn  State  University. 

The  measured  structural  transfer  functions  provide  the  relationship  between  point  forces  and 
velocities  at  positions  on  the  panel.  Several  methods  of  obtaining  the  structural  transfer  func¬ 
tions  were  conducted,  each  method  varying  in  the  type  of  input  force  and  sensing  transducer 
used.  The  variation  in  the  input  and  sensing  transducers  was  performed  in  order  to  achieve 
the  needed  accuracy  in  transducer  positioning  and  measurement  repeatability  and  to  provide  a 
common  phase  reference  between  the  transfer  functions  as  required  by  SRI. 

The  measurements  were  performed  on  a  34x34x0.125  inch  aluminum  panel  which  is  par¬ 
tially  covered  with  an  unconstrained  damping  layer  to  provide  realistic  damping.  The  panel 
is  mounted  to  a  wall  opening  of  a  transmission  loss  facility  and  the  mounting  approximates 
clamped  boundary  conditions,  Figure  37.  Figure  35  shows  a  photograph  of  the  laser  vibrom- 
eter  side  of  the  panel,  while  Figure  36  shows  a  photograph  of  the  hammer/shaker  side  of  the 
panel.  The  panel  is  discretized  into  nxn  patches  of  equal  size,  the  comers  of  which  are  the 
positions  where  the  measurements  are  taken  (sensing  positions)  and  control  forces  are  applied 
(excitation  positions).  It  should  be  noted  that  the  panel  has  symmetric  boundary  conditions  en¬ 
abling  the  principle  of  reciprocity  to  be  used,  i.e.,  the  interchanging  of  excitation  and  sensing 
positions  will  result  in  identical  transfer  functions. 


7.1.1  Experiment  1  -  Frequency  Domain  Analysis 

The  goal  of  the  initial  experimental  investigation  was  to  obtain  a  set  of  transfer  functions 
that  could  describe  completely  the  vibrational  motion  of  the  plate.  These  transfer  functions 
were  obtained  using  a  modal  impact  hammer  as  the  input  force  and  an  accelerometer  as  the 
dynamic  sensor.  Figure  38.  The  accelerometer  was  mounted  to  the  panel  at  a  set  of  specified 
sensing  positions  and  the  modal  impact  hammer  (manual)  excited  the  panel  at  a  set  of  specified 
excitation  positions.  The  signals  from  the  two  transducers  were  sent  to  an  HP  Analyzer,  where 
the  structural  transfer  function  was  computed  and  stored. 

This  procedure  was  deemed  unexceptable  for  two  reasons.  First,  the  repeatability  of  the  ex¬ 
citation  force  produced  by  modal  impact  hammer  was  questionable  due  to  human  error,  i.e. 
the  impact  was  excited  manually,  and  thus  the  excitation  force  was  not  always  normal  to  the 
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panel.  Second,  it  was  believed  that  the  mass  of  the  accelerometer  altered  the  panels  response 
when  mounted  to  the  panel  at  different  sensing  positions. 


7.1.2  Experiment  2  •  Frequency  Domain  to  Time  Domain  Analysis 

A  second  experimental  investigation  was  performed  to  eliminate  the  previously  mentioned 
problems  and  to  reduce  the  amount  of  recorded  data  by  allowing  the  input  force  to  remain 
at  a  single  position  while  only  moving  the  sensor  to  set  a  sensing  locations.  The  experimen¬ 
tal  transfer  functions  were  obtained  using  a  uni-morph  piezoelectric  shaker  as  the  force  input 
and  a  Polytech  laser  vibrometer  as  the  sensor,  Figure  39.  The  shaker  was  attached  to  a  force 
transducer,  which  was  in  turn  mounted  to  the  panel  at  a  specified  position.  The  HP  Analyzer 
provided  a  random  noise  signal  to  the  shaker.  This  method  of  excitation  allowed  for  a  repeat- 
able  input  force.  The  laser  vibrometer  is  a  non-contact  sensor  and  therefore  eliminates  the 
mass  loading  effects  of  the  accelerometer.  The  laser  vibrometer  was  manually  moved  to  set  of 
sensing  positions  in  order  to  measure  to  response  of  the  panel. 

It  was  determined  that  the  manual  positioning  of  the  laser  vibrometer  did  not  provide  the 
measurement  position  accuracy  required.  Consequently,  a  motorized  positioning  system  for 
the  laser  vibrometer  was  developed,  allowing  for  accurate  sensor  positioning,  and  the  data  was 
retaken.  The  use  of  the  mounted  shaker  and  motorized  positioning  of  the  laser  vibrometer  pro¬ 
duced  the  required  accuracy  in  measurement  repeatability  and  transducer  positioning  required 
in  developing  the  multichannel  control  circuit. 

It  was  later  determined  by  SRI  that  when  these  transfer  functions  were  converted  from  the  fre¬ 
quency  domain  to  the  time  domain,  there  was  a  loss  of  phase  reference  between  each  transfer 
function.  The  loss  of  phase  reference  is  inherent  when  using  random  excitation. 


7.1.3  Experiment  3  -  Frequency  Domain  to  Time  Domain  Analysis 

In  experiment  3,  an  electronic  impulse  hammer  was  used  as  the  excitation  force  because  it 
provided  both  a  common  phase  reference  between  transfer  functions  and  a  repeatable  excita¬ 
tion  force.  The  electronic  impulse  hammer  was  mounted  to  rigid  frame  located  adjacent  to 
the  panel.  Figure  40.  The  mounting  procedure  was  designed  such  that  it  was  easy  to  manually 
move  the  impulse  hammer  to  any  excitation  position.  The  laser  vibrometer  was  positioned  at 
a  set  of  sensing  positions  using  the  motorized  positioning  system. 

Since  this  experimental  set-up  achieved  the  needed  accuracy  in  transducer  positioning  and 
measurement  repeatability  and  provided  a  common  phase  reference  between  the  transfer  func¬ 
tions,  it  was  requested  by  SRI  that  data  was  taken  for  multiple  excitation  and  sensing  positions 
(similar  to  experiment  1,  except  a  complete  description  of  the  vibrational  motion  of  the  plate 
was  not  required). 
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Figure  35:  Laser  vibrometer  side  of  panel 


7.1.4  Conclusion 

A  means  of  obtaining  experimental  structural  transfer  functions  of  a  plate  was  developed 
in  order  to  provide  SRI  with  realistic  data  needed  to  design  a  multichannel  control  circuit. 
The  development  of  an  experimental  set-up  to  obtain  the  structural  transfer  functions  was  de¬ 
scribed.  The  experimental  set-up  met  the  required  accuracy  in  measurement  repeatability  and 
transducer  positioning  and  also  provided  a  common  phase  reference  between  the  individual 
transfer  functions. 


7.2  Validation  of  Green’s  Function  Properties  Using  PSU  DATA 

The  data  taken  by  PSU  was  compared  against  model  data  using  the  modal  plate  model.  Careful 
adjustment  of  the  model  parameters  yielded  a  reasonable  agreement  between  predicted  and 
measured  responses.  This  provide  us  confidence  in  using  the  data  to  perform  tests  in  validating 
our  control  approach.  Figure  41  shows  the  magnitude  and  phase  of  the  response  along  64 
evenly  spaced  points  along  a  strip  of  the  plate  due  to  point  forcing  at  a  single  location.  The 
response  as  a  function  of  frequency  is  shown  for  both  the  model  prediction  as  well  as  the 
measured  response  of  the  laser- vibrometer  output. 

In  order  to  validate  our  wavelet  transfer  function  compression  approach,  we  examined  the 
PSU  data  for  intrinsic  properties  of  the  response  data  to  impulsive  excitations.  In  particular, 
as  shown  in  [3],  the  impulse  response  of  a  structure  should  exhibit  discontinuities  at  certain 
derivatives  at  the  location  of  the  impulse.  This  behavior  enables  the  use  of  wavelets  to  com¬ 
press  the  transfer  function  from  inputs  to  outputs.  Figure  42  shows  the  velocity  response 
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Figure  36:  Hammer/shaker  side  of  panel 
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Figure  37:  Diagram  of  Panel  Mounted  in  Wall  Opening  of  Transmission  Loss  Facility 
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Figure  38:  Diagram  of  Experiment  1  Set-up 
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Figure  39:  Diagram  of  Experiment  2  Using  Laser  Vibrometer 
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Figure  42:  Velocity  response  due  to  force  input  using  model  data. 

along  a  strip  of  the  plate  (off  center)  due  to  a  force  input  at  a  single  location.  The  velocity 
at  64  evenly  spaced  points  are  shown  due  to  the  force  input,  using  model  data.  Note  that  the 
response  peaks  more  and  more  at  the  location  of  the  forcing  point  as  a  function  of  the  order 
of  the  derivative.  This  peaking  as  a  function  of  derivative  order  is  precisely  the  property  is 
exploited  by  wavelet  transforms  with  vanishing  moments  in  compressing  transfer  functions 
from  inputs  to  outputs. 

Figure  43  shows  the  result  of  measuring  the  velocity  at  these  64  evenly  spaced  points  due 
to  a  shaker  force  input.  The  top  plot  shows  the  raw  response  measured  by  PSU,  while  the 
subsequent  plot  shows  the  result  of  taking  derivatives  (using  finite  differences)  based  on  this 
raw  data.  This  plot  shows  excellent  validation  of  the  wavelet  transfer  function  compression 
approach,  since  the  peaking  property  is  clearly  exhibited  in  the  measured  data. 


7.3  Wavelet  Compression  Using  PSU  Data 

In  this  section  we  report  results  on  using  the  PSU  data  to  show  wavelet  compression  of  the 
neutralization  transfer  function.  In  this  case,  the  coupling  from  actuator  input  to  sensor  output 
along  a  strip  on  the  plate  was  used.  As  in  the  case  of  the  derivative  test,  64  evenly-spaced 
sensor  points  were  used.  Only  16  evenly-spaced  actuation  points  were  used,  due  to  limitations 
in  collecting  the  data.  These  16  were  chosen  to  be  collocated  with  the  first  16  sensor  points. 
Thus,  the  transfer  function  measured  has  16  inputs  and  64  outputs.  Figure  44  shows  the 
magnitude  of  the  transfer  function  at  a  single  frequency  and  the  magnitude  of  the  wavelet 
transform  at  that  frequency.  Note  that  the  “finger”-like  patterns  are  evident  in  the  wavelet 
transform  of  the  transfer  function.  Since  we  only  had  access  to  the  first  16  inputs,  rather  than 
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Figure  43:  Velocity  response  due  to  force  input  using  PSU  measured  data. 

the  full  64,  we  could  only  seek  to  verify  that  the  pattern  exhibited  matched  the  subset  of  the 
full  finger  pattern  corresponding  to  a  64  input  64  output  system. 

To  quantify  how  well  the  prediction  of  the  pattern  based  on  our  physical  smoothness  notions 
matched  the  measured  data,  we  applied  the  our  predicted  pattern  mask  to  the  measured  transfer 
function,  i.e.,  we  zeroed  out  the  coefficients  that  were  not  included  in  the  mask.  The  result  is 
shown  in  Figure  45. 

To  summarize  the  performance  of  the  compression  achieved  using  this  mask,  we  reconstructed 
the  transfer  function  in  physical  coordinates  and  determined  the  relative  error  to  be  only  .017 
percent.  The  original  and  reconstructed  transfer  function  can  be  visually  compared  in  Fig¬ 
ure  46. 
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Figure  44:  (1)  Transfer  function  magnitude  in  physical  coordinates  (r)  Magnitude  of  transfer  function 
in  wavelet  coordinates. 
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Figure  45 :  (1)  Predicted  mask  of  significant  coefficients  (r)  Mask  applied  to  measured  transfer  function. 
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Figure  46:  (1)  Transfer  function  magnitude  in  physical  coordinates  (r)  Reconstructed  transfer  function 
magnitude  in  physical  coordinates. 
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8  TECHNOLOGY  TRANSFER 


The  purpose  of  developing  control  technology  for  smart  materials  is  to  provide  enabling  tech¬ 
nology  that  will  be  applied  to  systems  with  large  numbers  of  sensors  and  actuators.  In  this 
section  we  summarize  past  and  continuing  efforts  to  effect  this  transfer  from  research  to  ap¬ 
plications. 

There  are  two  areas  of  development  in  which  efforts  at  technology  transfer  has  occurred;  trans¬ 
ducer  selection  and  transfer  function  compression.  As  discussed  in  Section  5,  the  transducer 
selection  or  placement  results  of  this  program  are  not  restricted  to  systems  with  large  num¬ 
bers  of  transducers,  so  that  work  has  been  tested  more  and  is  closer  to  application  that  the 
compression  results.  Among  the  applications  have  been:  a  Navy  structural  acoustic  control 
program,  the  SPICES  program,  and  a  NASA  aircraft  cabin  quieting  program.  The  application 
to  SPICES  is  the  primary  example  in  [1 1],  which  is  attached  as  Appendix  B. 

Since  the  “break-even”  point  for  the  compression  methodology  is  somewhere  on  the  order 
of  64  transducer  pairs,  a  considerably  greater  investment  in  transducers  is  required  before  an 
application  or  demonstration  can  benefit  from  using  our  results  instead  of  alternative  control 
technology.  Recent  and  current  demonstrations  and  programs  have  typically  focussed  on  de¬ 
velopment  of  individual  transducers  or  manufacturing  technology.  For  reasons  of  maximizing 
the  tangible  results  of  a  fixed  amount  of  funding  these  issues  have  received  priority.  However, 
technology  such  as  KIKO  will  be  necessary  to  take  these  efforts  beyond  “toy”  demonstra¬ 
tions,  and  it  is  imperative  that  KIKO  be  ready  when  the  applications  are  ready  for  that  step. 
We  summarize  the  additional  work  we  believe  needed  to  be  appropriately  ready  in  Section  9. 

It  appears  that  a  first  application  may  well  be  acoustic  noise  control.  We  have  had  numer¬ 
ous  discussions  with  groups  working  on  such  applications  in  aircraft  cabin  quieting  (the  same 
NASA  project  mentioned  above),  attenuation  of  acoustic  noise  for  rocket  payloads,  underwa¬ 
ter  acoustics  (also  in  the  DARPA  Smart  Materials  program),and  the  DARPA  MEMS  program 
(SRI’s  acoustic  noise  suppressor  program).  In  order  to  assist  in  control  development  (and  pos¬ 
sible  KIKO  appliation)  for  the  DARPA  Smart  Wing  program,  in  which  around  130  actuators 
are  anticipated  for  a  full  scale  system,  we  have  a  consulting  role  in  the  current  phase  of  that 
program.  We  also  believe  that  a  variation  of  this  technology  may  be  applicable  to  our  current 
work  on  the  ONR  Power  Electronic  Building  Block  (PEBB)  program.  (This  is  a  considerable 
departure  from  the  partial  differential  equation-based  theoretical  justification  in  [3]  included  in 
Appendix  C,  but  the  tremendous  proliferation  of  transducers  anticipated  in  the  use  of  PEBBs 
requires  some  form  of  transfer  function  compression  and  hierarchical  control,  and  the  wavelet 
approach  may  well  provide  the  right  set  of  coordinates  in  which  to  reduce  the  complexity  of 
models  for  the  desired  and  required  interconnections.) 

We  recognize  that  continuing  efforts  are  required  to  bring  this  work  to  its  useful  application, 
and  we  believe  that  KIKO  will  figure  prominently  into  future  programs  as  demonstrations  are 
developed  requiring  the  large  numbers  of  transducers  for  which  KIKO  is  intended. 
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9  CONCLUSIONS 


A  great  variety  of  military  and  commercial  applications  have  anticipated  requirements  for 
feedback  control  systems  using  large  arrays  of  sensors  and  actuators,  with  numbers  of  trans¬ 
ducers  on  the  order  of  100, 1000,  or  more.  These  applications  include  acoustic  noise  cancel¬ 
lation,  vibration  suppression,  adaptive  optics  and  control  of  turbulence  in  fluids. 

The  KIKO  methodology  comprises  a  procedure  for  designing  a  controller  that  is  locally  inter¬ 
connected  in  a  multiscale  hierarchy,  and  a  high-level  processor  architecture  for  implementing 
the  controller  in  a  scalable  manner.  A  key  feature  of  KIKO  is  that  it  provides  a  way  to  imple¬ 
ment  globally  optimal  control  designs  via  hierarchical  processing  and  interconnection.  This 
is  in  contrast  to  other  locally  interconnected  control  designs,  which  can  only  provide  local 
optimization  and  cannot  address  global  criteria. 

The  design  method  takes  as  input  transfer  functions  which  are  assumed  to  be  measured  on 
the  system  using  sensor  and  actuator  arrays.  Based  on  these  measured  transfer  functions  the 
method  uses  efficient  numerical  algorithms  to  design  a  controller  that  is  locally  interconnected 
in  a  banded  multiscale  representation.  The  high-level  processor  architecture  is  based  on  the 
fast  wavelet  transform  implemented  in  a  distributed  manner.  An  approach  for  integration  of 
multirate  techniques  into  the  multiscale  methodology  has  been  developed. 

The  design  procedure  has  been  tested  on  closed-loop  control  simulations  of  beams  and  plates, 
and  experimental  work  with  a  plate  confirmed  the  validity  of  key  assumptions.  Transducer 
placement  work  focussed  on  heuristic  subset  selection  from  arrays,  although  with  guaranteed 
performance  bounds. 

Although  long  term  application  requirements  suggest  major  extensions  of  KIKO  (say,  to  adap¬ 
tive  and  nonlinear  versions),  in  order  to  efficiently  exploit  and  extend  KIKO  developments  to 
date,  further  near-term  work  is  needed  to  ensure  that  KIKO  technology  is  ready  for  near  term 
demonstrations  and  applications  that  will  need  it. 

The  most  pressing  near-term  need  for  KIKO  is  to  provide  a  closed-loop  control  demonstration 
of  its  capability.  Otherwise,  further  needed  development  might  wait  until  a  demonstration 
is  planned  for  other  purposes,  and  it  is  discovered  by  trying  and  failing  that  other,  existing, 
control  methodology  is  not  adequate  for  the  job. 

The  most  productive  next  direction  for  continuing  research  is  likely  to  be  the  investigation 
of  second  generation  wavelets.  This  would  build  strength  in  three  important  issues:  treat¬ 
ing  irregularly  spaced  transducers,  treating  bounded  domains  of  transducer  distribution,  and 
reducing  the  break-even  point  at  which  KIKO  has  an  advantage  over  conventional  approaches. 
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A  Euler-Bernoulli  Green’s  Function 


In  this  section,  we  derive  an  analytic  expression  for  the  Green’s  function  of  a  forced  Euler- 
Bemoulli  beam  with  clamped  ends.  Figure  47  shows  a  schematic  of  the  system.  The  beam  is 
described  by  the  well  known  differential  equation  with  clamped  boundary  conditions 

EI^-p.Sui^w  =  FS(x-y), 

w  =  w'  =  ^  at  a;  =  0,  L,  (56) 


where  E  is  Young’s  modulus,  I  is  the  rotary  inertia  of  the  beam  cross  section,  S  is  the  beam’s 
cross-sectional  area,  ps  is  the  density  of  the  beam,  u  is  temporal  frequency  and  w  is  the 
displacement  of  the  beam’s  neutral  surface  as  a  function  of  x.  The  beam  is  of  length  L  and 
has  thickness  h.  For  simplicity,  we  define  the  following  nondimensional  parameters 


_  w  _  X 


u  = 


IpsSL^ 


w  F  = 


FL^ 


El  Elh 

In  terms  of  the  nondimensional  variables  above.  Equation  56  becomes 

d'^w  o 


-  y 

y=h- 


(57) 


urw  —  FS{x  —  y), 

r _ 


dx^ 

w  =  w'  =  0  at  X  =  0, 1. 


(58) 


In  the  following,  we  will  use  the  nondimensional  form  above  where  we  drop  the  overbars  for 
convenience. 


L 


Figure  47:  Schematic  for  the  Euler-Bernoulli  beam. 


The  derivation  of  the  Green’s  function  proceeds  as  follows.  The  beam  is  conceptually  “split” 
into  two  sections  on  the  left  and  right  of  the  forcing  point.  The  solutions  for  each  side  satisfy 
the  original  differential  equation,  a  boundary  condition  on  one  side,  and  a  set  of  matching 
conditions  joining  the  solutions  at  the  forcing  point.  The  resulting  solution  is  a  weak  solution 
of  the  original  differential  equation  in  that  it  is  twice  continuously  differentiable  on  x  €  [0, 1] 
with  a  jump  discontinuity  at  x  =  y  in  the  third  derivative.  The  advantage  of  this  weak  solution 
form  over  the  traditional  series  solution  in  terms  of  eigenvalues  and  eigenfunctions  of  the 
Sturm-Louville  problem  is  that  in  computations,  the  series  solution  must  be  truncated  at  some 
eigenvalue,  giving  rise  to  error  akin  to  the  Gibb’s  phenomenon.  The  weak  solution  has  no 
such  truncation  error.  Furthermore,  the  discontinuity  in  the  third  x-derivative  of  the  Green’s 
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Figure  48:  Splitting  of  the  beam  into  two  domains. 


function  intuitively  reflects  the  effect  of  forcing  in  a  way  that  the  infinite  sum  of  eigenfunctions 
to  form  this  same  step  does  not. 


We  start  by  splitting  the  beam  into  two  parts  about  the  forcing  point, 
This  gives  for  w 


{wi  0  <  X  <y 
wn  y  <  X  <1  ’ 


as  shown  in  Figure  48. 

(59) 


with 


d^wj 

dx^ 


—  U>^Wl 


wj  =  w'j  =  0 


d'^wji 

dx^ 


—  UpWjI 


Wii  =  w'jj  =  ^ 


=  0  0  <  a;  <  y 

at  X  =  0  and 

=  0  y  <  X  <  1 

at  X  =  1. 


(60) 

(61) 


In  order  for  the  split  system  to  satisfy  the  original  equation  for  w,  a  set  of  matching  conditions 
at  X  =  y  must  be  satisfied.  These  are  given  by 


wi  =  mi,  w\  =  w'jj,  w'l  =  w'Ij,  w'j  -  w'jj  =  F  at  x  =  y.  (62) 

The  appearance  of  the  jump  in  the  third  derivative  can  be  derived  by  evaluating  Equation  58 
on  the  left  and  the  right  of  y  and  matching  impulses.  The  resulting  split  system  couples  two 
boundary  value  problems. 


The  general  solution  for  wj  and  wn  is  given  by 

wi{x)  =  Cl  cos(\4^x)  +  C2  sin(-v/wx)  +  czCOsh{yJiJ)x  4-  C4sinh(v^x), 
wn{x)  =  C5  cos(v^x)  +  C6  sin(v^x)  +  C7COsh(Vwx)  +  C8sinh(-v/wx).  (63) 

The  two  boundary  conditions  on  each  end,  x  =  0, 1  and  the  four  matching  conditions  at  x  =  y 
give  eight  equations  in  the  eight  unknown  c^’s,  i  e  [1, 8].  The  system  for  the  Ci’s  is  given  by 


Transducer  Placement  for  Broadband  Active  Vibration  Control 
Using  a  Novel  Multidimensional  QR  Factorization 

Larry  P .  Heck  Koorosh  Naghshineh  Julia  A.  Olkin 

SRI  International  Western  Michigan  University  SRI  International 

Menlo  Park,  CA  94025  Kalamazoo,  MI  49008  Menlo  Park,  CA  94025 


Abstract 

This  paper  advances  the  state  of  the  art  in  the  selection  of  minimal  configurations  of  sensors  and 
actuators  for  active  vibration  control  with  smart  structures.  The  method  extends  previous  transducer 
selection  work  by  (1)  presenting  a  umfied  treatment  of  the  selection  and  placement  of  large  numbers 
of  both  sensors  and  actuators  in  a  smart  structure,  (2)  developing  computationaUy  efficient  techniques 
to  select  the  best  sensor-actuator  pairs  for  multiple  unknown  force  disturbcinces  exciting  the  structure, 
(3)  selecting  the  best  sensors  and  actuators  over  multiple  frequencies  ,  and  (4)  providing  bounds  on  the 
performance  of  the  transducer  selection  algorithms.  The  approach  is  based  on  a  novel,  multidimensional 
extension  of  the  Householder  QR  factorization  algorithm  appKed  to  the  frequency  response  matrices  that 
define  the  vibration  control  problem.  The  key  features  of  the  algorithm  are  its  very  low  computational 
complexity,  and  a  computable  bound  that  can  be  used  to  predict  whether  the  transducer  selection 
algonthm  will  yield  an  optimal  configuration  before  completing  the  search.  Optimal  configurations  will 
result  from  the  selection  method  when  the  bound  is  tight,  which  is  the  case  for  many  practical  vibration 
control  problems.  This  paper  presents  the  development  of  the  method,  as  weU  as  its  application  in  active 
vibration  control  of  a  plate. 


1  INTRODUCTION 


Active  vibration  control  represents  the  state  of  the  art  in  reducing  unwanted  vibrations  in  structures.  A 
current  focus  of  research  in  the  area  is  the  development  of  smart  materials  that  can  automatically  sense 
and  then  control  vibrations  induced  from  external  disturbances.  Smart  materials  for  vibration  control  ap¬ 
plications  typically  consist  of  embedded  sensors,  actuators,  communication  channels,  and  processors  that 
have  been  programmed  to  implement  the  prescribed  control  law  with  the  available  transducers  However 
current  control  algorithms  cannot  effectively  deal  with  large  system  configurations  (many  sensors  and  actua¬ 
tors  resulting  in  too  many  interconnections).  Thus,  subsets  of  transducers  are  often  employed  to  achieve  the 
desired  control.  Many  of  these  presume  the  existence  of  a  control  configuration  (i.e.,  layout  of  the  sensors  and 
actuators  and  the  interconnections  between  them)  or  use  ad-hoc  methods  to  specify  one,  and  focus  instead 
on  the  design  of  the  best  control  law.  This  is  in  spite  of  the  fact  that  the  complexity  and  performance  of  a 
control  system  are  largely  determined  by  the  underlying  control  configuration  (Nett  and  Spang,  1987). 

Recent  work  in  transducer  selection  algorithms  can  be  broadly  grouped  into  two  approaches.  The  first 
approach  assumes  that  permissible  transducer  locations  are  a  continuous  function  in  spatial  location  Ap¬ 
proaches  to  solve  this  problem  utilize  numerical  optimization  methods  where  the  transducer  locations  are 
included  as  optimization  parameters  (Chen  et  al.,  1991;  Skelton  and  DeLorenzo,  1983;  Wang  et  al.  1991)  In 
general,  these  methods  have  high  computational  complexity,  and,  as  a  result,  have  been  limited  to  small-^ale 
transducer  selection  probleins,  focusing  on  either  actuator  or  sensor  selection  (but  not  both).  Specifically 
the  methods  have  been  applied  to  problems  of  selecting  tens  of  transducers  at  a  single  frequency. 

The  second  transducer  selection  approach  treats  the  problem  where  only  a  fixed  (finite)  number  of 
transducer  locations  are  permissible.  Approaches  to  solve  this  discrete  selection  problem  have  focused  on  the 
use  of  heuristic  search  methods  to  treat  the  combinatorial  explosion  of  candidate  transducer  configurations. 


Work  in  this  area  includes  efforts  by  Snyder  and  Hansen  (1990).  In  their  paper,  the  authors  focus  on  the 
sensor  selection  problem,  and  present  an  efficient  least  squares  technique  to  evaluate  the  performance  of 
each  sensor  configuration.  However,  the  paper  does  not  address  the  problem  of  efficiently  searching  through 
all  combinations  and  instead  employs  exhaustive  search.  As  a  result,  the  approach  was  limited  to  sensor 
selection  only,  at  a  single  frequency.  Ruckman  and  Fuller  (1993)  extended  the  work  by  Snyder  and  Hansen 
by  focusing  on  selecting  actuators  (instead  of  sensors)  and  suggesting  the  use  of  the  general  body  of  subset 
selection  methods.  The  subset  selection  methods  include  sequential  search  methods  (e.g.,  forward  and 
backward  selection),  which  are  computationally  efficient  but  can  yield  unpredictably  poor  results.  Finally, 
work  by  Baek  and  Elliott  (1995)  employed  a  class  of  heuristic  search  methods  known  as  natural  algorithms 
for  the  discrete  transducer  selection  problem.  The  complexity  of  the  search  methods  (simulated  annealing 
and  genetic  algorithms)  is  much  more  efficient  than  exhaustive  search,  but  significantly  higher  than  the 
sequential  search  methods.  While  the  methods  provide  a  useful  alternative  to  the  sequential  search  methods 
(outperforming  them  in  some  cases),  the  computational  complexity  of  the  methods  restrict  their  use  to 
small-scale  transducer  selection  problems. 

This  paper  presents  a  new  discrete  transducer  selection  method  for  designing  active  vibration  control  con¬ 
figurations.  First,  the  method  extends  previous  transducer  selection  work  by  presenting  a  unified  treatment 
of  the  selection  emd  placement  of  large  numbers  of  both  sensors  and  actuators  in  a  smart  structure.  Second, 
the  new  method  has  a  very  low  computational  complexity  (lower  than  sequential  forward  selection)  making 
it  useful  for  very  large-scale  sensor  and  actuator  selection  problems  (thousands  of  sensor-actuator  pairs) . 
Third,  while  all  the  approaches  described  above  choose  transducers  assuming  knowledge  of  the  location  of 
the  excitation  forces  on  the  structure  (and  often  limit  the  excitation  to  a  point  force),  the  method  described 
in  this  paper  does  not  require  this  knowledge.  In  this  way,  the  best  transducer  locations  can  still  be  selected 
for  applications  where  the  locations  of  the  disturbing  forces  are  not  known  during  the  design  of  the  trans¬ 
ducer  configuration.  Fourth,  the  method  presented  here  chooses  the  best  actuators  over  multiple  frequencies. 
Finally,  the  paper  presents  a  computable  bound  that  can  be  used  to  predict  whether  the  transducer  selection 
algorithm  will  yield  an  optimal  configuration  before  completing  the  search.  Optimal  configurations  will  result 
from  the  selection  method  when  the  bound  is  tight,  which  is  the  case  for  many  practical  vibration  control 
problems. 

This  paper  is  organized  as  follows.  The  next  section  briefly  defines  the  vibration  control  problem  and 
specifies  measures  of  performance  that  the  control  law  uses.  Section  3  formally  defines  the  transducer  se¬ 
lection  problem  in  the  context  of  this  control  law  and  proceeds  to  describe  the  new  transducer  selection 
approach.  We  also  describe  the  complexity  of  the  various  selection  methods,  showing  that  the  multidimen¬ 
sional  QR  algorithm  is  computationally  more  efficient  than  the  other  leading  selection  algorithms  (including 
sequential  forward  selection).  Finally,  Section  4  presents  results  from  application  of  this  approach  to  a  model 
of  a  smart  plate  (i.e.,  with  embedded  sensors  and  actuators). 

2  VIBRATION  CONTROL  PROBLEM 


In  recent  years,  a  variety  of  approaches  have  been  developed  for  structural  vibration  control.  A  common 
requirement  of  these  approaches  is  knowledge  of  transfer  functions  relating  disturbances,  sensors,  and  actu¬ 
ators  to  each  other  and  to  the  areas  on  the  structure  that  are  to  be  controlled.  Knowledge  of  the  transfer 
functions  can  be  achieved  through  measurements  (i.e.,  on-line  system  identification),  or  through  the  use  of 
structural  models  (e.g.,  finite  element  models).  For  a  broad  class  of  structures,  the  transfer  functions  can 
be  described  by  a  set  of  linear  time  invariant  (LTI)  differential  equations.  The  state-space  form  is  a  popular 
compact  representation  that  writes  these  equations  in  first-order  form,  i.e., 

X  =  Jx  -I-  Kn  (1) 

y  =  L-x.  +  Mm  (2) 

where  J,  K,  L,  and  M  are  state  matrix  variables,  u  is  a  vector  formed  by  stacking  the  disturbances 
and  Na  actuator  control  inputs,  x  is  an  element  state  vector  of  first  and  second  derivative  terms,  and  y  is  a 


vector  of  sensor  outputs. 

Taking  the  Fourier  transform  of  the  state-space  system  in  Eq.  (2)  yields  an  equivalent  representation  in 
terms  of  the  frequency  response  matrix  or  plant  H(u): 

y(u;)  =  H{u;)U{u,)  ,  (3) 

where 

H{u)  =  L{j(jl  -  J)~^K  +M.  (4) 

The  frequency  response  matrix  can  be  divided  into  four  submatrices  E,  A,  S,  and  C: 

^  >  [  5(a;)  C{u)  J  ’  (5) 

where  E  6  is  the  complex  disturbance  transfer  matrix  that  relates  spatially  distributed  forces  (or 

pressures)  at  the  disturbance  source  to  spatially  distributed  displacements  (or  velocities,  accelerations)  in 
the  desired  vibration  control  zone  (  structural  locations  where  minimal  vibrations  are  desired);  A  is 

the  complex  actuator  transfer  matrix  relating  control  inputs  (e.g.,  voltages)  for  all  k  candidate  actuators  to 
displacements  in  the  desired  control  zone;  5  €  is  the  sensor  transfer  matrix  relating  spatially  distributed 
forces  at  the  disturbance  source  to  I  measured  normal  displacements  at  the  candidate  sensor  locations  on 
the  vibrating  structure;  and  C  G  is  the  feedback  coupling  between  the  actuators  and  the  sensors. 

Figure  1  shows  a  block  diagram  of  the  vibration  control  problem.  The  responses  from  the  disturbance  to 
the  sensors,  5,  are  sent  to  the  compensator  W.  The  outputs,  WS,  are  sent  through  the  actuators,  A,  to  the 
vibration  control  zone.  In  this  particular  control  implementation,  a  feedforward  neutralization  filter,  N,  is 
used  to  “neutralize”  the  feedback  coupling,  C,  between  the  actuators  and  sensors.  This  technique  is’ called 
(^-parameterization,  which  has  been  used  successfully  in  vibration  control  problems  to  simplify  control  design 
(Flamm  et  al.,  1995).  Thus,  AWS  is  an  estimate  of  the  vibrations  that  reach  the  control  regions  directly, 
denoted  by  E.  The  problem,  then,  is  to  find  a  causal  compensator,  W,  that  minimizes  the  expression 

H  II  ^  111  subject  to  ^  H  WS  ||i  <  c  ,  (6) 

W=Wt  U}=Wl 

where  w/  a,nd  w*  are  the  low  and  high  cutoff  frequencies  of  the  control  bandwidth,  respectively,  and  c  is  a 
scalar  specified  by  the  designer.  The  constraints  in  this  problem  are  to  ensure  that  the  controller  is  physically 
realizable,  i.e.,  that  the  energy  in  the  actuator  inputs,  WS,  does  not  saturate  the  actuators,  and  that  the 
compensator  is  causal. 

Methods  currently  exist  to  solve  the  control  problem  in  Eq.  (6).  An  efficient  time-domain  method  for 
vibration  and  noise  control  problems  that  uses  a  conjugate-gradient  minimization  approach  to  find  finite- 
impulse  response  (FIR)  controllers  is  described  by  Flamm  et  al.  (1995).  The  system,  developed  by  Olkin 
et  al.  (1994)  ,  has  been  used  to  design  multiple-input/multiple-output  controllers  for  broadband  noise  and 
vibration  control  with  large  spatial  extent. 

3  TRANSDUCER  SELECTION 

Referring  again  to  Figure  1,  the  choice  of  the  actuators  will  change  the  size  and  nature  of  the  S  and  A 
response  matrices.  Practical  limitations  on  the  number  and  placement  of  transducers  often  means 
that  vibrations  in  the  desired  control  zone  and  disturbance  forces  cannot  be  directly  measured,  nor  can  the 
forces  be  exactly  reproduced.  Satisfaction  of  the  control  objective  will  therefore  depend  on  the  ability  of  the 
transducers  to  estimate  and  then  manipulate  these  vibrations  at  the  desired  control  zones. 
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3.1  Mathematical  formulation 

For  the  transducer  selection  problem,  the  modified  control  objective  using  reduced  numbers  of  actuators 
(Heck  and  Naghshineh,  1994)  is  given  as 


min 

vv.Ha.n. 


E 

w=ai/ 


-  A  Ho  ly  n,  5  II2  subject  to  ^  ||  W  11,  5  II2  <  c  , 

U}—Wi 


(7) 


with  the  minimization  completed  over  two  additional  parameters  H^  e  (AT,,  <  m)  and  H,  € 

jN.y.n  ^  where  is  the  set  of  all  ixj  binary-valued  matrices.  The  parameters  are  called  actuator 
and  sensor  selectors,  respectively.  The  columns  of  !!„  and  rows  of  11,  are  unit  basis  vectors  corresponding  to 
actuator  and  sensor  locations  (e.g.,  unit  basis  vector  U3  =  [0  0  1  0  0]^  corresponds  to  the  third  transducer 
position).  An  actuator  selector  matrix  given  as 


Ha 


1  0  ■ 
0  0 
0  1 


would  select  the  first  and  third  actuator  responses  (columns)  of  A. 


(8) 


3.2  Key  issues 

The  problem  in  Eq.  (7)  is  difficult  to  solve  directly  because  the  quadratic  cost  function  is  non-differentiable. 
This  arises  because  of  the  discontinuous  nature  of  the  selector  matrices:  either  a  transducer  is  selected, 
or  it  is  not.  As  a  result,  the  problem  cannot  be  solved  through  direct  extensions  of  the  gradient-based 
optimization  methods  (such  as  the  conjugate  gradient  method),  and  instead  has  been  addressed  through 
direct  enumeration  (search)  methods.  The  complexity  of  the  enumeration  method  is  approximately  the 
complexity  of  a  complete  controller  design  and  performance  assessment  for  one  configuration  multiplied  by 
the  number  of  possible  transducer  configurations.  The  number  of  possible  configurations  grows  extremely 
rapidly  as  the  complexity  of  the  system  increases. 

For  the  transducer  selection  problem,  given  m  candidate  sensors  and  n  candidate  actuators,  the  number 
of  distinct  subsystems  with  k  sensors  and  /  actuators  will  be 

rn\  n\ 

\  k  )  \  I  )  k\{m  —  k)\ l\{n  — 1)\ 

The  total  number  of  distinct  subsystems  will  be 

kzzl  1=1  ^  ^  ^  / 

where  each  subsystem  corresponds  to  a  selected  subset  of  the  candidate  transducers. 

Even  for  modest-size  problems,  the  number  of  candidate  transducer  configurations  in  Eq.  (10)  becomes 
exceedingly  large.  For  example,  the  number  of  possible  configurations  or  subsystems  in  a  system  with  40 
sensors  and  40  actuators  is  approximately  l.le-f23  (assuming  only  subsystems  with  equal  numbers  of  sensors 
and  actuators  are  considered).  To  put  this  into  perspective,  if  each  configuration  could  be  evaluated  in  1 
second  (i.e.,  a  complete  control  design  with  the  compensator  computed),  then  it  would  take  over  3.4e+13 
centuries  to  evaluate  every  possible  configuration!  As  a  result  of  the  large  numbers  of  possible  transducer 
configurations,  a  practical  method  for  transducer  selection  will  necessarily  consist  of  suboptimal,  but  highly 
efficient  measures  to  assess  controller  quality  along  with  fast  enumeration  algorithms. 
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3.3  Transducer  selection  approach:  single  frequency 

The  approach  to  transducer  selection  presented  in  this  section  employs  computationally  efficient  heuristic 
methods  that  avoid  the  high  combinatorics  described  above.  The  approach  consists  of  two  main  steps. 
First,  the  problem  shown  in  Eq.  (7)  is  broken  down  into  smaller,  more  tractable  problems.  Then,  after 
the  large  set  of  candidate  transducer  configurations  is  reduced  to  a  small  number,  full  controller  designs  are 
completed.  In  the  first  step,  the  approach  seeks  submatrices  of  A  and  S  selected  by  11^  and  TLs  that  extract 
the  non-redundant  components  of  the  original  matrices.  In  addition,  it  is  assumed  that  the  residual  error 
of  the  controller  with  these  full  matrices  can  be  made  small,  i.e.,  E  ^  AWS  at  all  frequencies.  Beginning 
with  the  single  frequency  case,  the  approach  approximates  the  controller  with  selected  transducers  such  that 
AUaWU^S  ^  AWS.  To  achieve  this,  a  subspace  distance  measure,  Dg,  can  be  employed  to  quantify  the 
difference  between  the  original  and  candidate  submatrices.  A  popular  measure  in  the  statistics,  numerical 
analysis,  and  signal  processing  literature  (Golub  and  Van  Loan,  1989)  is  given  as 

Ds{A,AUa)  =  sin0(7^(A),7^(A^a))  (11) 

D,(S,U,S)  =  sin0(7^(5’’),7^((^,5f)),  (12) 

where  Tl  denotes  range  space,  and  sin  0  is  the  sine  of  the  largest  principal  angle  between  the  range  subspaces 
of  the  transfer  function  matrices.  The  largest  principal  angle  between  two  subspaces  F  and  G  is  the  angle 
between  the  unit  vector  in  F  whose  orthogonal  projection  onto  G  gives  the  largest  residual.  The  sine  of 
the  largest  principal  angle  is  equivalent  to  the  2-norm  of  the  difference  (residual)  between  the  orthogonal 
projections  onto  the  range  spaces  of  the  matrices. 

Since  the  largest  principal  angle  corresponds  to  a  worst-case  residual,  it  can  be  argued  that  the  selector 
matrices  should  be  chosen  to  make  the  distance  measures  in  Eqs.  (11)  and  (12)  as  small  as  possible.  However, 
for  the  vibration  control  problem,  such  an  approach  could  lead  to  large  actuator  inputs,  which  would  violate 
the  constraints  in  Eq.  (7).  To  see  this,  consider  the  following  example  (Golub  and  Van  Loan,  1989).  Let 
5,  n,  =  /  in  Eq.  (7),  where  I  is  the  identity  matrix.  Also,  let  a;  =  so  that  the  minimization  is  completed 
over  a  single  frequency.  We  can  rewrite  the  minimization  problem  as 

min  \\E-ATlaW  \\l  subject  to  ||  W"  ||i  <  c  .  (13) 

WfJJa 

If  the  goal  is  to  find  the  two  best  actuators,  and 
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then  choosing  the  first  two  actuators  (first  two  columns  of  A)  gives  the  residual  \\  E  -  AUa  Hi  =  0?  but 
II  {AUaYE  II2  =  II  Wls  Hi  =  where  is  the  Moore-Penrose  pseudoinverse  operator  (Stewart, 

1984),  As  e  gets  smaller,  the  actuator  weights  will  grow  and  eventually  violate  the  constraint  in  Eq.  (13). 
On  the  other  hand,  the  other  two  possible  subsets  of  actuators  will  lead  to  small  weight  norms,  but  with  a 
much  worse  residual. 

This  example  highlights  the  tradeoff  that  exists  between  choosing  an  independent  set  of  actuators  and 
choosing  a  set  that  minimizes  the  residual  control  error.  A  similar  tradeoff  occurs  in  the  use  of  subset  selection 
methods  (Golub  et  al.,  1976)  to  solve  general  rank-deficient  least-squares  problems.  A  solution  developed 
for  this  problem  uses  a  Householder  QR  with  column-pivoting  factorization.  For  the  transducer  selection 
problem,  the  factorization  can  be  used  to  efficiently  find  an  independent  set  of  actuators  and  sensors.  In 
addition,  if  the  proper  number  of  transducers  is  selected,  the  residuals  in  Eq.  (7)  can  be  made  small. 

The  Householder  QR  algorithm  can  be  described  as  follows.  Let  G  E  matrix  be  representative  of 

the  A  and  S  matrices  in  Eqs.  (11)  and  (12).  The  QR  decomposition  of  G  can  be  written  as 


G 


where  Q  G  has  orthonormal  columns  =  /),  Rn  is  a  ^  x  k  matrix  with  a  condition  number 

approximately  equal  to  (Ti/cTk^  ||  R22  ||  is  of  the  order  (Tk+i^  and  k  is  the  rank  of  G,  Pivoting  (or  permuting) 
the  columns  of  the  QR  factorization  can  be  used  to  make  R22  small.  For  a  rank-deficient  matrix  with 
r  =  rank(G),  column  pivoting  can  yield 

Q-^ap  =  [  ]  ,  (15) 

where  iln  G  is  upper  triangular  and  nonsingular,  and  R12  €  With  the  column  partitionings 

GP  =  [flTpj ,  iTpj , . . . ,  and  <5  =  [91 , 92 ,  •  •  •  >  9m] ,  and  r,j  denoting  the  {i,  j)-th  element  of  R,  we  have 

min{r,k} 

9ph  -  ^  ''ikqi  €  span{9i,  92,  •  •  • ,  9r}  for  A;  =  1  :  n  .  (16) 

i=l 

This  implies  that  the  range(G)  =  span{gi,  ♦ .  .,?r}  =  span{^pi,  •  •  -  other  words,  column  pivoting 

can  be  used  to  find  a  subset  of  columns  in  a  rank-deficient  matrix  that  has  the  same  range  space  as  the 
original  matrix,  making  the  principal  angles  between  the  original  and  submatrix  go  to  zero.  In  addition,  the 
Householder  procedure  described  in  Golub  and  Van  Loan  (1989)  ensures  that  the  resulting  submatrix  of  G  is 
also  well-conditioned  (i.e.,  that  the  columns  are  sufficiently  independent).  Selecting  actuators,  for  example, 
with  this  procedure  (columns  of  A  in  Eq.  (7))  will  make  the  distance  measure  in  Eq.  (11)  go  to  zero  and 
yield  a  subset  with  the  same  control  authority  as  the  original,  densely  packed  set  of  actuators. 

In  most  vibration  control  problems  of  interest,  the  rank  of  the  A  and  matrices  will  not  be  exactly 
zero.  Rather,  because  of  numerical  approximations,  the  rank  can  only  be  specified  in  terms  of  a  threshold 
on  singular  values,  called  the  numerical  rank  (Stewart,  1984).  A  consequence  is  that  nonzero  differences  will 
exist  between  the  range  spaces  of  the  original  and  selected  submatrices.  As  a  result,  a  method  is  required  to 
determine  the  circumstances  under  which  the  differences  between  the  range  spaces  is  small.  A  useful  bound 
(Golub  and  Van  Loan,  1989)  can  be  determined  that  specifies  these  circumstances,  i.e., 

sin0(yl,  Allfl)  <  (Tr+i  II  II  (17) 

sin 6(5^,  (Sn,)^)  <  (T,+i  II  ii^iMl  .  (18) 

where  and  o-g+i  are  the  (r  +  l)st  and  (s-f  l)st  largest  singular  values  of  A  and  S',  respectively,  and  the 
submatrix  Ru  (this  in  general  will  be  a  different  matrix  for  A  and  S)  in  the  QR  factorizations  is  guaranteed 
to  be  of  the  order  for  A  and  for  S.  These  bounds  show  a  sufficient  condition  for  small  errors;  i.e., 
if  there  exists  a  large  gap  in  successive  singular  values,  called  the  spectral  gap,  then  the  error  will  be  small. 
For  transducer  selection,  these  bounds  provide  a  useful  criterion  for  the  optimality  of  a  given  configuration. 
Also,  they  give  a  direct  correspondence  between  the  numerical  rank  of  the  frequency  response  matrices  and 
the  number  of  required  transducers;  e.g.,  if  the  numerical  rank  of  the  actuator  frequency  response  matrix  is 
A:,  then  k  actuators  are  required. 

SA  Broadband  extension 

The  goal  for  the  broadband  transducer  selection  problem  is  to  minimize  the  expression  in  Eq.  (7)  with  a 
single  configuration  over  all  frequencies  of  interest.  For  the  broadband  case,  we  extend  the  Householder 
procedure  to  compute  factorizations  over  multiple  matrices.  The  approach  is  as  follows.  Assume  for  some  k 
and  frequency  LJi  that  we  have  computed  Householder  matrices  /fi, . . . ,  and  permutations  Pi, . . . , 
such  that 

(Hk^iioJi)  •  -  Hi{u;i))G{uJi){Pi{u;i)  •  •  •Pfe_i(a;i)) 
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where  G  and  R^22~^\u)i)  €  The 

single-frequency  Householder  method  moves  the  column  of  G  corresponding  to  the  largest  column  in  R22{^i) 
to  the  A:th  column  position  (the  lead  position  in  the  i^22(^i)  matrix),  where  largest  is  specified  by  the 
selection  of  a  vector  norm.  The  column  exchanges  are  accomplished  with  the  permutation  matrix 
This  is  followed  by  a  left  multiplication  by  a  Householder  matrix  that  zeroes  all  the  subdiagonal  components 
of  the  largest  R22{^i)  column  vector.  The  procedure  stops  when  Ar  —  1  =  rank(G(wi)). 

For  the  broadband  case  (multiple  G  matrices),  instead  of  permuting  a  G(u;i)  based  on  the  column  with 
the  largest  vector  norm  in  permute  G(u;i), . , . ,  G(u;n)  based  on  one  of  eight  possible  performance 

measures  on  the  matrix  consisting  of  a  column  from  each  ii22(^i),  •  •  • ,  R22{^n)  matrix. 

These  performance  measures,  discussed  and  analyzed  in  more  detail  in  Olkin,  Heck,  and  Naghshineh 
(1996)  ,  can  be  selected  depending  on  the  application.  The  first  four  performance  measures  are  the  three 
induced  matrix  norms,  /i,  I2  and  loo  and  the  Frobenius  norm.  The  fifth  performance  measure  chooses  the 
smallest  singular  value  of  each  actuator  which  helps  minimize  the  effect  of  deep  nulls  in  an  actuator  frequency 
response,  which  can  lead  to  nonrobust  control  designs  with  poor  performance.  The  sixth  measure  chooses 
the  actuator  response  with  the  largest  average  value  over  all  frequencies  and  response  locations.  The  seventh 
measure  is  the  ratio  of  the  minimum  to  maximum  actuator  response  over  frequency,  reflecting  the  fact  that 
flatter  actuator  frequency  responses  are  often  desirable  due  to  noise  control  robustness  and  performance 
considerations.  The  last  performance  measure  attempts  to  minimize  the  effect  of  spatial  nulls  by  choosing 
the  maximum  row  sum  of  the  actuator  matrix. 

As  was  the  case  with  a  single  matrix,  this  permutation  step  is  followed  by  left  multiplying  each  G(a;i), 
...,  G{oJn)  with  the  corresponding  Householder  matrices  that  will  zero  the  subdiagonal  components  of  the 
selected  columns  for  each  i^22(wi), « *^,R22{^n)  matrix.  Figure  2  illustrates  this  multidimensional  broadband 
extension  of  the  Householder  QR  factorization  with  column  pivoting.  Note  that  P  is  not  a  function  of 
frequency;  hence  this  matrix  of  ones  and  zeros  is  the  same  for  each  frequency. 

For  the  broadband  case,  the  error  bounds  in  Eqs.  (17)  and  (18)  hold  for  each  frequency,  i.e., 

sin0(A(wi),i4(wi)na)  <  o-r+i(wi)  II  II 

sin0(A(w2),>l(w2)na)  <  0-r+l(w2)  II  iill(W2)“Ml 

sin0(^(w„),^(fa;„)na)  <  (rr+i(w„)  ||  || 

and 

sin0(5(a;i)^,(5(wi)n^)^)  <  o-,+i(a)i)  ||  i?ii(wi)“^  || 

sin0(5(w2)^,  (5(a;2)n,)^)  <  cr,+i(«2)  ||  ^11(^2)“^  || 

sin0(5(c4)„)^,(5(w„)n,)^)  <  <r,+i(w„)  ||  ||  . 

The  bounds  hold  for  each  frequency  because  the  derivation  of  the  bounds  does  not  depend  upon  the  specific 
properties  of  the  QR  factorization.  Rather,  the  proof  holds  for  a  general  Ru  matrix.  However,  the  quality 
of  the  bound  may  be  adversely  affected  if  ||  R'[^  ||  is  not  on  the  order  of  <7^^. 

3.5  Computational  complexity  of  approach 

The  algorithm  for  implementing  the  transducer  selection  approach  described  above  is  computationally  effi¬ 
cient  and,  as  a  result,  can  be  used  for  very  large-scale  problems  involving  hundreds  to  thousands  of  sensors 
and  actuators.  Specifically,  suppose  we  are  choosing  k  actuators  out  of  Nd  total  actuators,  with  Nq  error 
sensors,  over  Nf  frequencies.  Each  of  the  following  transducer  selection  methods  requires  some  sort  of  com¬ 
putation  of  the  residual  of  the  active  vibration  control  system.  Let  G (residual) a?  represent  the  flop  count 
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for  computing  the  system  residual  with  k  actuators.  Then  Table  1  gives  a  general  complexity  overview  of 
some  actuator  selection  methods,  and  the  dominating  term  by  way  of  easy  comparison.  The  methods  are  (i) 
Exhaustive  Search,  in  which  all  combinations  of  k  actuators  are  considered,  (ii)  Best  A:,  in  which  each  actua¬ 
tor  is  evaluated  individually  and  the  Best  k  are  chosen,  (iii)  sequential  forward  selection,  in  which  actuators 
are  added  to  the  chosen  set  based  on  which  new  actuator  will  behave  the  best  with  those  already  chosen, 
and  finally,  (iv)  Multidimensional  QR,  the  algorithm  presented  in  this  paper.  Choosing  k  out  of  Ng  sensors 
yields  similar  relative  complexity  counts  as  in  Table  1,  with  Multidimensional  QR  still  more  computational 
efficient  than  the  SFS  method.  For  brevity,  we  include  results  for  the  actuators  only. 


Table  1:  Algorithm  Complexity  For  Choosing  k  out  of  Nd  Actuators 


Method 

General  Complexity 

Dominating  Term 

Exhaustive  Search 

1v^3=(6A^:FT20F) 

Multidim.  QR 

i^^(column-selection)-f 

AT/C?  (Householder  transform) 

P{NfN^  +  N^)  +  kNfNf  +  NjNe 

Seq.  Forward  Sel. 

ELi  -  *  +  l)<!?(residual)j} 

k^NfNdNe  +  k'^NfNd  -  k^Nj 

Best  k 

NdNf  0  (residual)  i 

It  can  be  difficult  to  compare  the  dominating  terms,  so  as  way  of  example,  consider  a  system  with  40 
possible  actuators  in  which  only  20  can  be  kept.  Assume  there  are  20  error  sensors  measuring  the  actuator 
responses  over  300  frequencies.  Then  the  complexity  for  these  four  methods  are  shown  in  Table  2,  printed 
in  order  of  optimality  of  solution. 


Table  2:  Example  Complexity  Comparison 


Method 

Complexity 

Exhaustive  Search 

6.9el9 

Sequential  Forward  Selection 

3.8e9 

Multidimensional  QR 

9.9e7 

Best  k 

1.7e6 

Notice  that  the  Multidimensional  QR  algorithm  is  more  efficient  than  sequential  forward  selection,  and 
comparable  with  Best  K.  In  addition,  the  method  is  guaranteed  to  perform  as  well  or  better  than  both  SFS 
and  Best  K.  As  will  be  demonstrated  in  the  next  section,  the  efficiency  of  this  procedure  enables  practical 
transducer  selection  for  large-scale  problems. 


4  NUMERICAL  RESULTS 

The  methods  described  in  the  previous  section  were  applied  to  the  problem  of  controlling  vibrations  on  a 
smart  plate  model  (embedded  sensors  and  actuators).  The  SPICES  consortium’s  solid-plate  model  (Flamm 
et  al.  1995)  was  considered. 

As  shown  in  Figure  3,  the  plate  is  a  layered  composite  material  46  x  46  x  1.2  cm.  (18x18x0. 5  in.), 
connected  to  an  infinite,  rigid  base  on  the  bottom  of  the  plate  at  the  four  corners  (indicated  by  the  black  dots). 
The  plate  was  excited  with  a  point  force  with  x,  y,  and  z  components.  The  force  was  applied  to  the  plate 
through  a  stiff  tripod  (all  three  connections  with  the  plate  move  together  to  simulate  connectors  between 
the  plate  and  attached  machinery).  The  forcing  function  was  an  impulse  with  fiat  frequency  response.  The 
first  268  natural  frequencies  and  mode  shapes  of  this  3858-node  plate  model  were  used  in  the  computation 
of  frequency  response  matrices  (admittance)  relating^^citation  force  to  displacement. 


The  sensor  models  used  in  this  numerical  experiment  measure  point  displacements.  For  the  sensor 
placement  studies,  361  evenly  spaced  candidate  sensor  locations  were  used.  The  candidate  locations  were  on 
a  19  X  19  grid  (2.5  cm=l  in.  from  nearest  neighbor)  in  the  middle  layer  (denoted  by  the  dashed  line).  For  the 
actuator  placement  studies,  simple  two-dimensional  actuators  were  modeled  and  incorporated  into  the  plate 
model.  The  actuators  have  length  and  depth  (into  the  plate),  but  no  width.  Figure  4  shows  a  cross-section 
side  view  of  the  plate  with  an  embedded  actuator.  The  directional  forces  in  the  model  are  shown  as  ±F,  which 
produce  an  effective  bending  moment  about  the  center  of  the  actuator.  A  total  of  180  candidate  actuator 
locations  were  specified  for  the  actuator  selection  studies.  From  a  top  view  of  the  plate  (xy  projection), 
the  actuator  locations  formed  a  checkerboard  pattern,  as  illustrated  in  the  right  plot  in  Figure  3.  The  lines 
between  the  dots  on  the  checkerboard  pattern  represent  the  actuators,  and  all  horizontal  and  vertical  lines, 
represented  by  letters  in  this  schematic,  are  considered  actuators.  The  goal  of  the  vibration  control  system 
was  to  reduce  the  transmission  of  forces  over  500-4000  Hz  from  the  tripod  excitation  to  the  connectors  at 
the  four  corners  of  the  plate.  Thus,  the  four  corners  of  the  plate  model  represent  the  vibration  control  zone. 
All  control  law  designs  used  FIR  implementations  designed  using  the  conjugate-gradient  method  described 
in  Section  2.  The  FIR  weights  were  found  by  solving  the  minimization  problem  in  Eq.  (6). 

To  select  the  sensors  and  actuators,  the  continuous,  time-domain,  state-space  model  responses  were 
sampled  and  used  to  compute  the  frequency  response  matrices  S  and  A  in  Eq.  (7).  The  sensor  and  actuator 
selction  was  then  performed  using  these  matrices.  Finally,  the  FIR  controller  was  implemented  and  the 
controlled  system  responses  computed. 

The  S  matrix  consisted  of  361  rows  (corresponding  to  sensor  responses),  and  the  A  matrix  had  180 
columns  (actuator  responses) .  The  broadband  extension  of  the  Householder  QR  with  column-pivoting  fac¬ 
torization  was  applied  to  these  matrices  to  find  the  columns  of  A  and  rows  of  S  that  formed  matrices  with 
approximately  the  same  range  space  as  the  original  matrices.  As  described  in  Section  3.4,  the  matrix  2-norm 
was  used  by  the  Householder  algorithm  to  select  the  columns  of  the  response  matrices.  The  numerical  rank 
of  the  S  matrix  (with  a  threshold  set  @  -30  dB)  was  approximately  3.  Based  on  this,  three  sensors  were 
chosen.  Their  locations  correspond  to  the  connector  points  between  the  tripod  and  the  plate  in  Figure  3. 
The  numerical  rank  of  the  A  matrix  was  approximately  4.  The  left  plot  in  Figure  5  shows  the  locations 
of  the  actuators  and  sensors  chosen  by  the  broadband  QR  method.  The  right  plot  shows  the  uncontrolled 
(dashed)  and  controlled  (solid)  time-domain  impulse  responses  for  the  control  bandwidth  at  the  x  =  17  in., 
y  =  17  in.  connector.  As  can  be  seen,  this  control  configuration  quickly  reduces  the  vibration  transmission 
from  the  tripod  to  the  connectors.  Figure  6  shows  the  power  spectral  density  (psd)  of  the  uncontrolled  and 
controlled  vibration  responses  for  the  sensor-actuator  layout  chosen  by  the  multidimensional  QR  algorithm. 
The  lines  shown  represent  the  root  mean  squared  (RMS)  of  the  values  of  the  responses  at  all  four  feet  from 
the  X,  y,  and  z  tripod  excitation.  As  can  be  seen,  the  controller  effectively  reduces  the  transmission  power 
over  the  control  bandwidth  (500-4000  Hz), 

To  compare  the  QR-based  designs  with  manually  selected  actuator  locations,  additional  designs  were 
completed,  chosen  to  be  intuitively  competitive.  The  results  are  shown  in  Figures  7  through  10  (all  cases 
used  the  QR  sensor  configuration).  As  can  be  seen,  the  performance  for  each  design  were  inferior  to  those 
of  the  QR-based  actuator  configuration.  Table  3  shows  the  in-band  control  performance  (500-4000Hz)  at 
each  of  the  four  feet  (or  corners)  of  the  plate  in  dB.  The  three  configurations  of  actuator  placement  are  the 
QR-based  configuration  and  the  two  comparison  configurations  shown  in  Figures  7  and  9.  Notice  that  the 
performance  at  all  four  corners  for  the  QR-based  configuration  is  about  the  same  level,  whereas  there  is 
much  more  variation  for  the  other  two  configurations,  and  even  some  enhancement. 


5  SUMMARY 


This  paper  has  presented  a  computationally  efficient,  unifying  method  for  broadband  sensor  and  actuator 
selection.  The  method  is  based  on  the  Householder  QR  with  column-pivoting  factorization  of  the  system’s 
frequency  response  matrices.  Using  the  Householder  algorithm,  upper  bounds  were  presented  on  the  errors 
induced  by  selecting  a  subset  of  transducers.  The  bounds  can  be  computed  efficiently,  as  they  are  in  terms 
of  the  singular  values  of  the  response  matrices.  Als^^  they  serve  as  sufficient  conditions  for  optimality  of 


Table  3:  Performance  Comparison. 


Configuration 

In-Band  Performance  (dB) 

Corner  1 

Corner  2 

Corner  3 

Corner  4 

QR 

10.9 

10.9 

11.4 

11.4 

Figure  7 

6.86 

-.24 

.76 

-.40 

Figure  9 

10.6 

5.8 

2.3 

3.2 

the  transducer  selection  method,  such  that  rank-deficient  response  matrices  with  large  spectral  gaps  in  the 
singular  values  cause  the  bounds  in  Eq.  (17)  and  Eq.  (18)  to  reach  equality  instead  of  inequality.  Finally,  the 
multidimensional  extension  of  the  Householder  QR  procedure  developed  in  this  paper  provides  a  practical 
method  of  selecting  sensors  and  actuators  that  yield  effective  control  over  a  broad  band  of  frequencies.  The 
new  method  is  more  efficient  than  sequential  forward  selection  and  has  an  order  of  complexity  approaching 
that  of  simply  choosing  the  best  K  transducers  independently.  However,  the  new  method  performs  as  well 
or  better  than  both  competing  algorithms.  Future  work  will  extend  this  method  to  minimize  the  required 
interconnections  between  the  sensors  and  actuators. 


6  ACKNOW^LEDGEMENTS 

This  work  was  supported  by  ARPA  Contracts  MDA972-93- 2-0010  and  N00014-94-C-0176.  The  authors  wish 
to  thank  the  SPICES  team  for  providing  the  plate  model  used  in  the  numerical  experiments  presented  here, 
and  David  Flamm,  Bill  Nowlin,  Paul  Titterton,  and  Ken  Chou  for  helpful  discussion  related  to  this  work. 

7  REFERENCES 

K.  Baek,  S.  Elliott,  1995,  “Natural  algorithms  for  choosing  source  locations  in  active  control  systems,” 
ASME  J.  Sound  and  Vib.^  vol.  186,  no.  2,  pp.  245-267. 

G.  Chen,  R.  Bruno,  and  M.  Salama,  1991,  “Optimal  placement  of  active/passive  members  in  truss 
structures  using  simulated  annealing,”  AIAA  Journal^  vol.  29,  no.  8. 

R.  Clark  and  C.  Fuller,  1992,  “Optimal  placement  of  piezoelectric  actuators  and  polyvinylidence  fluoride 
error  sensors  in  active  structural  acoustic  control  approaches,”  J.  Acoust.  Soc.  Am.,  vol.  92(3),  pp.  1521- 
1533. 

D.  Flamm,  L.  Heck,  P.  Titterton,  W.  Nowlin,  J,  Olkin,  and  K.  Chou,  1995,  “Control  system  design 
for  the  spices  smart  structure  demonstrations,”  in  SPIE  Conference  on  Smart  Structures  and  Materials, 
pp.  237-248. 

G.  H.  Golub,  V.  Klema,  and  G.  W.  Stewart,  1976,  “Rank  degeneracy  and  least  squares  problems,”  Tech. 
Rep.  TR-456,  Depart,  of  Computer  Science,  University  of  Maryland,  College  Park,  MD. 

G.  H.  Golub  and  C.  F.  Van  Loan,  1989,  Matrix  Computations.  Baltimore:  The  Johns  Hopkins  University 
Press,  second  ed, 

L.  Heck  and  K,  Naghshineh,  1994,  “Large-scale,  broadband  actuator  selection  for  active  noise  control,” 
in  Proceedings  of  Nat.  Conf  in  Noise  Control  Engineering,  (Ft.  Lauderdale,  FL). 

C.  Nett  and  H.  Spang,  1987,  “Control  structure  design:  A  missing  link  in  the  evolution  of  modern  control 
theories,”  in  American  Control  Conference,  vol.  29. 

J.  A.  Olkin,  M.  S.  Freed,  and  P.  D.  Jungers,  1994,  “SRI  weights  algorithm  and  performance  simulation 
(SWAPS)  code.,”  tech,  rep.,  SRI  International. 

J.  Olkin,  L.  Heck,  and  K.  Naghshineh,  1996,  “Automated  placement  of  transducers  for  active  noise 
control:  Performance  measures,”  in  International  Conference  on  Acoustics,  Speech,  and  Signal  Processing, 
(Atlanta,  GA).  C.  Ruckman  and  C.  Fuller,  1993  “Optimizing  actuator  locations  in  feedforward  active  control 

86 


systems  using  subset  selection,”  in  Supplement  to  2nd  Conf.  on  Recent  Advances  in  Active  Control  of  Sound 
and  Vibration,  pp.  S122-S133. 

K.  Skelton,  R.E.  DeLorenzo,  M.L.,  1983,  “Selection  of  noisy  actuators  and  sensors  in  linear  stochastic 
systems,”  Journal  of  Large  Scale  Systems,  Theory,  and  Applications,  vol.  4,  pp.  109-136. 

S.  D.  Snyder  and  C.  H.  Hansen,  1990,  “Using  multiple  regression  to  optimize  active  noise  control  system 
design,”  in  J.  of  Sound  and  Vib.. 

G.  W.  Stewart,  1984,  “Rank  degeneracy,”  SIAM  J.  Sci.  Statist.  Comput,  vol.  5,  pp.  403-413. 

B.  Wang,  R.  Burdisso,  and  C.  Fuller,  1991,  “Optimal  placement  of  piezoelectric  actuators  for  active 
control  of  sound  radiation  from  elastic  plates,”  Proceedings  of  Nat.  Conf  in  Noise  Control  Enaineerina 
pp.  267-274. 

B.  Wang,  1996,  “Optimal  placement  of  microphones  and  piezoelectric  transducer  actuators  for  far-field 
sound  radiation  control,”  J.  Acoust.  Soc.  Am.,  vol.  99,  no.  5. 


87 


Table  Legends 


Table  1 
Table  2 
Table  3 


:  Algorithm  Complexity  For  Choosing  k  out  of  Nd  Actuators. 
:  Example  Complexity  Comparison. 

:  Performance  Comparison. 
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Figure  Legends 


Figure  1:  Block  diagram  of  the  vibration  control  problem.  The  responses  from  the  disturbance  to  the 
sensors,  S',  are  sent  to  the  compensator  W.  The  outputs,  WS,  are  sent  through  the  actuators,  A,  to 
the  vibration  control  zone.  A  feedforward  neutralization  filter,  A,  is  used  to  “neutralize”  the  feedback 
coupling,  (7,  between  the  actuators  and  sensors. 

Figure  2:  Multidimensional  Householder  QR  factorization  with  column  pivoting. 

Figure  3:  The  left  plot  shows  the  solid-plate  model  developed  by  the  SPICES  consortium.  The  right  plot 
shows  the  pattern  for  the  180  candidate  actuator  locations  embedded  in  the  SPICES  plate. 

Figure  4:  Actuator  Model. 

Figure  5:  Locations  of  the  actuators  and  sensors  chosen  by  the  broadband  QR  method  (left).  The  uncon¬ 
trolled  (dashed)  and  the  controlled  (solid)  time-domain  impulse  responses  for  the  control  bandwidth 
for  the  connector  located  at  x  =  17zn.,  t/  =  17m.  (right). 

Figure  6:  Plot  of  the  RMS  values  of  the  uncontrolled  (solid)  and  the  controlled  (dashed)  responses  at  all 
four  feet  for  the  configuration  chosen  by  the  QR  algorithm. 

Figure  7:  Locations  of  the  actuators  and  sensors  chosen  by  the  authors  for  comparison  to  the  QR  method 
(left).  The  uncontrolled  (dashed)  and  the  controlled  (solid)  time-domain  impulse  responses  for  the 
control  bandwidth  for  the  connector  located  at  x  =  17m.,  y  =  17m.  (right). 

Figure  8:  Plot  of  the  RMS  of  the  uncontrolled  (solid)  and  the  controlled  (dashed)  responses  at  all  four  feet 
for  the  configuration  chosen  by  the  authors. 

Figure  9:  Locations  of  the  actuators  and  sensors  chosen  by  the  authors  for  comparison  to  the  QR  method 
(left).  The  uncontrolled  (dashed)  and  the  controlled  (solid)  time-domain  impulse  responses  for  the 
control  bandwidth  for  the  connector  located  at  x  =  17m.,  y  =  17m.  (right). 

Figure  10:  Plot  of  the  RMS  of  the  uncontrolled  (solid)  and  the  controlled  (dashed)  responses  at  all  four 
feet  for  the  configuration  chosen  by  the  authors. 
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Figure  1:  Block  diagram  of  the  vibration  control  problem.  The  responses  from  the  disturbance  to  the  sensors, 
S,  are  sent  to  the  compensator  W .  The  outputs,  WS,  are  sent  through  the  actuators.  A,  to  the  vibration 
control  zone.  A  feedforward  neutralization  filter,  N,  is  used  to  “neutralize^’  the  feedback  coupling,  C,  between 
the  actuators  and  sensors. 


Z  (inch«s) 


Figure  5:  Locations  of  the  actuators  and  sensors  chosen  by  the  broadband  QR  method  (left).  The  uncontrolled 
(dashed)  and  the  controlled  (solid)  time-domain  impulse  responses  for  the  control  bandwidth  for  the  connector 
located  at  x  —  17m,,  y  —  YJin.  (right). 
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Figure  6:  Plot  of  the  largest  singular  values  of  the  uncontrolled  (solid)  and  the  controlled  (dashed)  responses 
at  all  four  feet  for  the  configuration  chosen  by  the  QR  algorithm. 
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Figure  7:  Locations  of  the  actuators  and  sensors  chosen  by  the  authors  for  comparison  to  the  QR  method 
(left).  The  uncontrolled  (dashed)  and  the  controlled  (solid)  time-domain  impulse  responses  for  the  control 
bandwidth  for  the  connector  located  at  x  =  17m.,  y  =  17m.  (right). 
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Figure  8:  Plot  of  the  largest  singular  values  of  the  uncontrolled  (solid)  and  the  controlled  (dashed)  responses 
at  all  four  feet  for  the  configuration  chosen  by  the  authors. 
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Figure  9;  Locations  of  the  actuators  and  sensors  chosen  by  the  authors  for  comparison  to  the  QR  method 
(left).  The  uncontrolled  (dashed)  and  the  controlled  (solid)  time-domain  impulse  responses  for  the  control 
bandwidth  for  the  connector  located  at  x  —  17m.,  y  =  17m.  (right). 
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Figure  10:  Plot  of  the  largest  singular  values  of  the  uncontrolled  (solid)  and  the  controlled  (dashed)  responses 
at  all  four  feet  for  the  configuration  chosen  by  the  authors. 
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Abstract 

In  this  paper  we  analyze  the  representation  of  integral  operators  whose  kernels  are 
Green’s  functions  for  a  class  of  linear  differential  equations  using  wavelets  with  a  finite 
number  of  vanishing  moments.  In  particular,  we  show  how  wavelets  can  be  used  to  generate 
a  sparse  representation  of  these  operators.  We  show  that  the  matrix  associated  with  the 
discretized  integral  operator  represented  in  the  wavelet  basis  is  sparse  and,  in  particular, 
contains  multiple  bands  of  various  widths.  The  particular  banded  structure  of  the  wavelet 
representation  of  the  operator  follows  from  the  fact  that  the  associated  Green’s  function 
is  smooth  away  from  the  source  point  and  is  singular  at  some  order;  i.e.,  for  some  T,  its 
Tth  derivative  is  discontinuous  at  the  source  point.  We  derive  bounds  on  the  magnitude 
of  the  coefficients  of  the  integral  operator  in  the  wavelet  basis  as  a  function  of  scale  and 
position  and,  in  particular,  in  terms  of  whether  or  not  the  coefficient  lies  within  a  band. 
Based  on  these  bounds,  we  can  approximate  the  operator  by  ignoring  coefficients  not  lying 
within  these  bands,  thus  producing  a  sparse  representation.  This  sparse  representation  is 
extremely  beneficial  for  numerical  applications  in  which  one  would  like  to  apply  the  Green’s 
function  operator  efficiently:  normally  if  such  an  operator  mapped  N  points  into  N  points, 
it  would  require  0{N^)  operations;  however,  with  the  wavelet  transform,  the  mapping 
would  require  only  0((4M  -1-  2jLM  -f  27(1  —  j)M  —  S)N)  operations,  where  2M  is  the 
length  associated  with  the  support  of  the  wavelet  function,  L  =  log2  N  —  1,  and  7  = 

An  application  example  in  which  this  is  important  is  the  control  of  Smart  Structures  in 
which  a  large  number  of  embedded  sensors  and  actuators  must  be  coordinated  in  order  to 
achieve  disturbance  rejection  on  the  surface  of  the  structure. 

Keywords:  wavelets,  Green’s  functions,  control,  structures 
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1  Introduction 


The  literature  on  the  use  of  wavelets  to  efficiently  represent  signals  and  mathematical  operators 
has  become  extensive,  spanning  a  wide  range  of  applications.  While  much  of  the  mathematical 
aspects  of  wavelet  representations  of  operators  has  its  origins  in  the  work  of  Meyer  [9],  the  work 
of  Beylkin  et  al.  [1]  was  perhaps  the  first  presentation  of  numerical  results  that  generated  a 
great  deal  of  interest  in  the  idea  of  using  wavelets  to  develop  sparse  representations  of  operators. 
In  [1],  wavelets  were  used  to  compress  Calderon-Zygmund  and  pseudo-differential  operators.  In 
this  work  it  was  shown  that  representations  of  these  operators  using  tensor  products  of  one¬ 
dimensional  wavelets  were  extremely  efficient  in  the  sense  that  few  wavelet  coefficients  were 
needed  to  represent  the  operators  with  arbitrarily  specified  accuracy.  In  this  paper,  we  focus  on 
the  use  of  wavelet  transforms  in  compressing  integral  operators  whose  kernel  is  a  Green’s  function 
associated  with  a  linear  differential  equation.  While  there  has  been  much  work  in  using  wavelets 
to  solve  partial  differential  equations  (PDEs)  numerically  using  finite  difference  and  Galerkin 
methods  [3,  5],  our  objective  is  to  present  a  novel  analysis  of  Green’s  functions  associated  with 
the  integral  form  of  the  solution  to  a  linear  differential  equation  and  to  show  how  the  sparsity  of 
wavelet  representations  of  Green’s  function  integral  operators  follows  naturally  from  properties 
of  both  wavelets  and  Green’s  functions. 

Our  paper  presents  an  analysis  of  the  use  of  wavelets  with  vanishing  moments  to  represent  inte¬ 
gral  operators  associated  with  Green’s  functions  for  linear  differential  equations.  We  show  that 
integral  operators  of  this  type  can  be  efficiently  represented  in  such  wavelet  bases,  i.e.  very  few 
wavelet  coefficients  need  to  be  retained  to  faithfully  represent  the  operator.  Furthermore,  there 
exist  extremely  efficient  algorithms  for  computing  the  wavelet  transform  [4],  making  the  overall 
transformation  scheme  a  viable  method  for  efficiently  realizing  the  Green’s  function  operator. 
This  efficiency  is  important  for  such  applications  as  numerical  integration,  the  numerical  solu¬ 
tion  of  differential  equations,  the  solution  of  inverse  problems,  as  well  as  for  ongoing  work  we  are 
performing  in  the  area  of  the  control  of  “Smart  Structures”  [2].  Furthermore,  we  are  currently 
working  on  extensions  of  this  work  to  the  solution  of  PDEs. 

We  apply  our  analysis  to  the  class  of  Green’s  functions  associated  with  differential  equations  of 
the  following  form: 

Cu  =  f  (1) 

where  £  is  a  linear,  differential  operator,  u  is  the  solution  to  the  differential  equation,  and  /  is 
the  forcing  term.  The  solution  u  can  be  represented  in  integral  form  as  follows. 

u(r)  =  f  G{r,s)f{s)ds  (2) 

j  s 

where  G(r,  s)  is  the  Green’s  function  kernel  of  the  integral  operator.  We  will  refer  to  the  vaxiable 
r  as  the  observation  variable,  and  the  variable  s  as  the  source  variable. 

The  main  contribution  of  this  paper  is  an  analysis  of  the  structure  of  the  representation  of  Green’s 
function  integral  operators  in  wavelet  bases.  In  particular,  we  characterize  the  sparsity  of  the 
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transformed  operator  in  terms  of  basic  properties  of  Green’s  functions.  Our  analysis  relies  on 
the  following  fundamental  property  of  Green’s  functions:  For  linear  differential  operators  with 
smooth  coefficients,  the  associated  Green’s  function  is  smooth  (in  the  C°°  sense)  in  both  the  source 
and  observation  variables,  except  at  points  where  they  are  equal.  Furthermore,  at  points  where 
they  are  equal,  i.e.,  where  the  source  and  observation  locations  coincide,  for  each  coordinate  in 
both  observation  and  source  variables,  the  highest  derivative  of  the  Green’s  function  with  respect 
to  that  coordinate  is  singular.  This  fact  reflects  the  basic  physical  property  that  a  Dirac  delta 
forcing  function  introduces  a  singularity  in  the  highest  order  derivative  of  the  Green’s  function 
with  respect  to  the  observation  variable  at  the  point  of  forcing. 

Our  paper  presents  an  analysis  of  the  structure  of  the  Green’s  function  integral  operator  coeffi¬ 
cients  in  the  wavelet  representation,  and  in  particular,  this  analysis  reveals  explicitly  the  sparse 
structure  of  the  operator  in  this  representation  as  a  function  of  spatial  scale  and  position,  and 
as  a  function  of  the  structure  of  the  singularities  in  the  Green’s  function.  Our  analysis  exploits 
the  fact  that  wavelets  are  particularly  well-suited  to  analyzing  singularities.  The  use  of  wavelet 
transforms  in  analyzing  singularities  is  widely  recognized  in  such  applications  as  image  analysis 
[7],  where  wavelets  help  reveal  the  structure  of  edges  in  an  image.  In  [7],  wavelets  were  used 
to  characterize  the  Lipschitz  regularity  of  a  function  based  on  the  decay  rate  of  the  wavelet 
coefficients  as  a  function  of  scale.  In  [8]  wavelets  axe  used  to  analyze  singularities  arising  in 
turbulence.  In  all  these  applications,  the  fact  that  wavelets  have  a  finite  number  of  vanishing 
moments  is  essential,  enabling  the  wavelets  to  discriminate  between  smooth  portions  of  functions 
and  their  points  of  singularity.  Note  that  the  thrust  of  this  paper  is  to  provide  a  theoretical  basis 
for  approximating  Green’s  function  integral  operators  as  sparse  operators  in  the  wavelet  domain. 
An  example  application  where  this  is  particularly  useful,  which  we  will  not  discuss  in  detail  in 
this  paper,  is  the  control  of  Smart  Structures.  In  this  application,  one  might,  for  example,  want 
to  control  the  vibration  of  the  structure  using  embedded  sensors  and  actuators.  In  [2],  a  method 
is  used  in  which  the  transfer  function  from  actuator  inputs  to  sensor  outputs  is  modeled  as  a 
Green’s  function  operator  and  an  extremely  efficient  algorithm  is  developed  to  control  a  structure 
containing  large  arrays  of  sensors  and  actuators. 

The  paper  is  organized  as  follows.  In  Section  2,  we  give  background  on  orthonormal  wavelets 
with  compact  support.  Note  that  while  we  perform  our  analysis  specifically  using  orthonormal 
wavelets,  our  results  do  not  in  fact  rely  on  this  property,  as  we  will  point  out  later.  We  then 
review  the  vanishing  moments  property  of  wavelets  and  some  results  on  analyzing  the  singularities 
of  functions  in  terms  of  Lipschitz  regularity.  In  Section  3,  we  develop  our  characterization  of 
the  wavelet  representation  of  Green’s  function  integral  operators,  using  fundamental  properties 
of  Green’s  functions  as  well  as  the  vanishing  moments  property  of  wavelets.  We  also  provide 
bounds  for  the  magnitude  of  the  coefficients  as  a  function  of  scale  and  position  of  the  wavelet 
basis  functions  associated  with  each  coefficient.  Then,  in  Section  4  we  give  numerical  examples 
of  our  results  based  on  a  differential  equation  describing  the  behavior  of  a  flexible  beam  (Euler- 
Bernoulli  beam).  These  examples  illustrate  the  multiband  structure  of  the  discretized  integral 
operator  in  the  wavelet  basis  as  well  as  the  scale  and  position  dependencies  of  the  operator 
coefficients  predicted  by  our  results.  Finally,  in  Section  5  we  state  conclusions  of  our  work  and 
point  out  extensions  that  we  are  currently  researching. 
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2  Background 


In  this  section  we  give  background  on  wavelets  with  an  arbitrary  number  of  vanishing  moments. 
We  treat  specifically  the  case  of  orthonormal  wavelets  of  compact  support,  although  our  results 
can  be  developed  for  the  case  of  non-orthogonal  wavelets  as  well.  The  issue  of  compact  support 
of  the  basis  functions  is,  however,  a  key  issue  in  our  development.  We  then  give  some  background 
on  the  relationship  between  the  number  of  vanishing  moments  of  a  wavelet  basis  function  and 
the  ability  of  the  basis  set  to  represent  smooth  functions  efficiently. 


2.1  Orthonormal  Wavelets  of  Compact  Support 

The  basic  structure  of  scale-based  wavelet  transforms^  builds  on  the  idea  that  a  space  of  functions, 
such  as  can  be  spanned  by  a  basis  set  consisting  of  dilations  and  translations  of  a  single, 

properly  chosen  function: 

—  n)  for  m,n  integer}  (3) 

The  basic  wavelet  function  ^  is  in  general  oscillatory  and  compressed  in  a:  -  hence  the  name. 
For  our  purposes  one  of  the  most  interesting  aspects  of  bases  of  the  type  in  eq.(3)  is  the  fact 
that  under  certain  conditions,  transforming  a  signal  in  these  basis  functions  can  be  interpreted 
as  performing  what  is  known  as  a  multiresolution  analysis  of  the  signal  [6] .  In  this  analysis  the 
signal  is  represented  at  different  resolutions  using  a  scaling  function,  $,  which  is  used  to  blur  the 
signal  at  different  resolutions. 

The  multiresolution  analysis  interpretation  of  the  wavelet  transform  can  be  summarized  by  the 
following  key  points. 


•  Translations  of  the  scaling  function  for  a  fixed  scale  m,  {$(2"* a;  —  n)  for  n  integer},  define 

a  subspace  of  At  this  scale,  a  signal  is  represented  by  projecting  the  signal  onto 

this  subspace  . 

•  Translations  of  the  wavelet  function  for  a  fixed  scale  m,  {^(2™ a:  —  n)  for  n  integer},  span 
the  subspace  defined  by  the  difference  between  the  scaling  subspace  at  scale  m  and  the 
scaling  subspace  at  scale  m  —  1. 

•  Wavelet  basis  coefficients  can  be  computed  recursively  in  scale,  resulting  in  an  extremely 
efficient  procedure. 


In  [4]  a  class  of  wavelets  was  introduced  in  which  the  wavelets  are  compactly  supported  and 
orthonormal.  Furthermore,  a  characteristic  of  this  class  of  wavelets  is  that  their  transforms  can 

^As  opposed,  for  example,  to  time-frequency  representations  such  as  Gabor  transforms,  which  some  include  in 
the  wavelet  class. 
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^  c(m-2,n) ,  d(m-2,n) 


c(m-1,n),d{rT>-1,n) 
c(m,n) ,  d(m,n) 


Figure  1:  Domain  of  Lattice  Algorithm  for  Computing  Wavelet  Coefficients 


be  computed  efficiently  via  a  scale  recursion.  Let  us  define  the  scaling  coefficients  and  wavelet 
coefficients  as  follows  for  some  signal  f{x): 

/OO 

$(2”‘x  —  n)f(x)dx  (4) 

-OO 

/OO 

—  n)f{x)dx  (5) 

-OO 

It  is  shown  in  [4]  that  the  following  equations  can  be  used  to  compute  the  scaling  and  wavelet 
coefficients  recursively  in  scale  from  fine  to  coaxse: 


c(m,n) 

=  y^c(m  +  l,k)h(2n 

k 

-k) 

(6) 

d{m,  n) 

=  '^c{m  +  l,k)g{2n 

-k) 

(7) 

k 


Each  of  the  recursions  in  eq.’s(6,7)  can  be  interpreted  as  discrete  filtering  followed  by  downsam¬ 
pling  by  a  factor  of  two.  The  scaling  function  filter  h{n)  is  related  to  the  wavelet  filter  g{n)  by: 
g{n)  =  (— l)”/i(l  —  n).  Furthermore,  for  compactly  supported  wavelets  h{n)  and  g(n)  are  FIR 
filters  (typically,  of  length  <  12).  The  computational  complexity  for  the  algorithm  described  by 
eq.’s(6,7)  is  order  0{NM)  where  2M  is  the  length  of  h  and  g  and  N  is  the  number  of  coefficients 
at  the  initial  (finest)  scale.  The  algorithm  flow  is  illustrated  in  Figure  1  for  the  case  of  4-tap  FIR 
filters  h  and  g,  where  the  c  and  d  coefficients  are  computed  recursively  up  the  lattice. 

The  inverse  transform,  i.e.  reconstructing  the  scaling  coefficients  at  the  finest  scale  from  wavelet 
coefficients  at  all  coarser  scales,  consists  of  the  following  recursion  in  scale  from  coarse  to  fine, 
which  represents  the  adjoint  of  the  operation  described  in  eq.’s(6,7). 

c(m  -1- 1, 2n)  =  ^  c(m,  k)h{n  —  k)  +  ^  d{m,  k)g{n  —  k)  (8) 

k  k 

The  algorithm  flow  would  correspond  to  propagating  down  the  lattice  in  Figure  1.  Again,  the 
computational  complexity  is  0{NM). 

What  these  fast  wavelet  algorithms  suggest  is  the  possibility  of  using  fast  basis  transformations 
in  applying  our  Green’s  function  integral  operator,  and,  assuming  this  operator  is  sparse  in  the 
wavelet  basis,  for  implementing  a  scheme  which  would,  overall,  be  extremely  efficient.  Note  that 
although  we  are  dealing  with  an  operator  that  can  be  viewed  as  a  superposition  integral,  it  is 
in  general  not  shift-invariant;  thus,  the  discretized  operator  would  not  be  implementable  using 
an  FFT.  Also,  the  wavelet  transform  is  highly  parallelizable,  offering  the  possibility  of  extremely 
fast  computations  using  multiprocessor  schemes. 
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2.2  Compression  Properties  of  Wavelets 


In  this  section  we  review  some  results  in  [1]  and  [7]  on  using  wavelets  both  to  characterize  the 
regularity  of  functions  and  to  compress  or  “sparsify”  operators.  These  results  motivate  our  use 
of  wavelet  bases  to  compress  operators  of  the  form  in  eq.(2).  We  also  derive  useful  results  on 
Lipschitz  regularity  of  a  function  at  a  point  and  its  relationship  to  polynomial  approximation  of 
the  function  at  that  point.  These  results  will  aid  in  our  analysis  of  representing  Green’s  function 
integral  operators  in  wavelet  bases. 

The  main  idea  in  using  wavelets  to  compress  operators  is  to  take  advantage  of  the  fact  that 
the  wavelet  bases,  for  example  those  designed  in  [4],  can  be  used  to  represent  smooth  operators 
efficiently.  The  feature  of  these  wavelets  that  allows  such  compression  is  referred  to  as  the 
property  of  “vanishing  moments,”  which  is  described  as  follows. 

/OO 

'^'(a:)a:’"da:  =  0for  m  =  0,l,...M-l  (9) 

-OO 


This  property  is  a  direct  result  of  the  fact  that  the  Fourier  transform  of  the  filter  h{n)  has  M 
zeros  at  a;  =  TT  and  that,  consequently,  the  Fourier  transform  of  the  filter  g{n)  has  M  zeros  at 
a;  =  0.  It  is  also  true  that  the  number  of  vanishing  moments,  M,  is  proportional  to  the  support 
length  of  the  wavelet  or  the  length  of  the  discrete  filters  h{n)  and  g{n). 

The  reason  this  property  is  essential  to  compressing  smooth  operators  is  the  fact  that  when  one 
integrates  a  smooth  function,  f{x),  against  a  wavelet  basis  function,  'F(2’”a;— n),  satisfying  eq.(9), 
the  first  M  terms  of  the  Taylor  expansion  of  /  vanish  and  one  is  left  with  the  integral  of  the 
product  of  ^  and  the  Mth  order  remainder  term  of  the  expansion.  What  this  means  for  integral 
operators  is  that  we  can  relate  compression  efficiency  of  the  wavelet  basis  in  terms  of  smoothness 
properties  of  the  kernel  of  the  operator.  In  [1]  the  class  of  Calderon-Zygmund  integral  operators 
are  considered,  which  have  the  property  of  being  smooth  away  from  the  diagonal.  In  particular, 
kernels  K{x,y)  associated  with  operators  in  this  class  satisfy  the  following: 

|8fif(x,!,)|+|a"iir(x,y)|  <  (11) 

From  the  fact  that  the  Mth  partials  of  K  decay  strongly  away  from  the  diagonals,  along  with 
the  vanishing  moment  property  in  eq.(9),  one  can  show  that  the  wavelet  coefficients  decay  in  a 
similar  fashion  as  a  function  of  the  distance  between  coefficients,  i.e.  the  distance  between  the 
centers  of  the  basis  functions  corresponding  to  two  coefficients. 


In  this  paper  we  consider  operators  that  are  not  restricted  to  the  Calderon-Zygmund  class.  In 
fact,  the  operators  can  be  frequency  dependent,  spatially  varying  in  their  parameters,  and  in 
general  depend  on  boundary  conditions.  Furthermore,  what  we  exploit  in  our  kernels  is  not 
the  fact  that  they  decay  away  from  their  diagonals;  in  fact  for  many  applications  they  do  not. 
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Rather,  we  exploit  the  fact  that  Green’s  function  kernels  exhibit  singularities  on  the  diagonals 
at  a  derivative  of  some  order. 

In  order  to  understand  more  precisely  the  implications  of  applying  wavelets  with  vanishing  mo¬ 
ments  to  a  kernel  with  a  singularity  on  its  diagonal,  we  state  precisely  what  we  mean  by  a 
singularity  of  a  function  at  a  point  at  some  order  derivative.  Furthermore,  we  will  be  exploiting 
the  vanishing  moments  property  of  wavelets  by  taking  Taylor  expansions  of  the  Green’s  func¬ 
tion  up  to  some  order.  Thus,  we  would  like  to  understand  the  relationship  between  the  order 
derivative  at  the  point  at  which  there  is  a  discontinuity  and  the  order  of  the  Taylor  expansion 
about  that  point.  In  the  following,  we  derive  statements  about  polynomial  approximations  of  a 
function  about  a  point  at  which  we  know  there  is  a  singularity  at  some  order  derivative. 

Let  us  define  the  local  smoothness  of  a  function  at  a  point  in  terms  of  Lipschitz  exponents  in  the 
following  way  [7]: 


Definition  2.1  For  n  a  positive  integer,  a  function  f{x)  is  Lipschitz  n  +  I  at  xq  if  there  exist 
constants  A  and  ho  >  0,  and  a  polynomial  of  order  n,  Pn{x),  such  that  for  0  <  h  <  ho 

\f{xo  +  h)-P„{h)\<Ah^+'^  (12) 


We  now  state  the  following  theorem  (see  Appendix  A  for  a  proof),  which  will  allow  us  to  show 
precisely  why  the  limited  differentiability  of  a  function  /  prevents  us  from  approximating  /  by 
a  polynomial  to  arbitrary  order. 


Theorem  2.1  For  n  a  positive  integer,  if  f{x)  is  Lipschitz  n  -)-  1  at  xo  and  f{x)  is  (7"  for 
0  <  |a:  —  iCol  <  7,  7  >  0,  then  f{x)  is  n  times  differentiable  at  xo  and  the  polynomial  Pn(h)  is 
equal  to  the  first  n  -f  1  terms  of  the  Taylor  series  expansion  of  f{x)  about  xo. 


Note  that  this  is  a  strong  theorem  relating  differentiability  of  a  function  and  approximation  of  a 
function  by  a  polynomial.  In  particular,  for  example,  if  the  nth  derivative  of  /  does  not  exist  at 
Xo,  there  is  no  polynomial  of  order  n  that  approximates  /  with  order  accuracy.  We  state 
this  more  precisely  in  the  following  corollaries,  the  first  of  which  is  essentially  a  restatement  of 
Theorem  2.1  in  contrapositive  form. 


Corollary  2.1  Let  f{x)  be  C°°  in  a  neighborhood  about  Xo,  excluding  Xo,  i.e.  C°°  for  x  such 
that  0  <  I®  —  xol  <  7;  7  >  0.  For  n  a  positive  integer,  if  f{x)  is  not  n  times  differentiable  at  Xo, 
then  f{x)  is  not  Lipschitz  n  -|-  1  at  xo.  Furthermore,  f  is  not  Lipschitz  j  at  rco  for  j  >  n  -f-  1. 
In  particular,  for  all  constants  A  and  ho  >  0,  and  for  all  polynomials  of  order  j,  Pj{x),  j  >  n, 
there  exists  an  h  such  that  0  <  h  <  ho  and 

\f{xo  +  h)-Pffh)\>AV+^  (13) 
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This  corollary  can  be  strengthed  into  the  following  corollary. 


Corollary  2.2  Let  f{x)  be  C°°  in  a  neighborhood  about  xq,  excluding  xq,  i.e.  C°°  for  x  such 
that  0  <  |a;  —  xo|  <  7,  7  >  0.  For  n  a  positive  integer,  if  f{x)  is  not  n  times  differentiable  at 
Xo,  then  for  all  constants  A  and  ho  >  0,  and  for  all  polynomials  of  order  j,  Pj{x),  j  >  n,  there 
exists  an  h  such  that  0  <  h  <  ho  and 

\f{xo  +  h)-Pj{h)\>Ah^+'^  (14) 


Proof 

We  know  from  Corollary  2.1  the  following.  For  all  constants  A  and  ho  >  0,  and  for  all  polynomials 
of  order  j,  Pj{x),  j  >  n,  there  exists  an  h  such  that  0  <  h  <  ho  and 

\f{xo  +  h)-Pj{h)\>Ah^+^  (15) 

Suppose  there  exists  &  j  >n  and  constants  B,ho>  0  such  that 

\f{xo  +  h)-Pj{h)\<Bh^+^  (16) 

for  0  <  ^  <  ho.  Let  Pj{h)  =  J2i-o  Pkh^  and  let  P„(/i)  denote  the  first  n  + 1  terms  of  Pffh).  Then, 
from  the  triangle  inequality  eq.(16)  can  be  written  as  the  following. 

\f{Xo  +  h)-Pn{h)\<  ^  +  (17) 

A:=n+1 


Then,  we  can  find  (7,  Aq  >  0  such  that 

\f{xo  +  h)-Pn{h)\<Ch'^^^  (18) 

for  0  <  <  ho.  But  this  contradicts  eq.(15)  for  the  case  j  =  n  and  Pj{h)  =  Pn{h). 

□ 

As  we  will  see,  these  corollaries  will  be  used  in  showing  how  limited  differentiability  of  the  kernel 
of  our  integral  operator,  i.e.  the  Green’s  function,  affects  our  ability  to  exploit  the  vanishing 
moments  property,  eq.(9),  of  wavelets. 


3  Representation  of  Green’s  Function  Operator  in  Wavelet 
Basis 


In  this  section  we  use  properties  of  Green’s  functions  to  characterize  the  structure  of  the  associ¬ 
ated  integral  operator  in  the  wavelet  basis.  This  transformed  operator  is  sparse  in  the  sense  that 
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the  magnitude  of  the  coefficients  along  various  bands  is  large  relative  to  the  magnitude  of  the 
coefficients  off  these  bands,  allowing  one  to  approximate  the  operator  by  ignoring  the  off-band 
coefficients  with  relatively  little  error.  We  show  explicitly  how  the  limited  differentiability  of  the 
Green’s  functions  along  their  diagonal  accounts  for  the  structure  of  these  bands.  Furthermore,  we 
give  bounds  on  the  relative  magnitude  of  the  coefficients  on  and  off  these  bands.  Note  that  while 
our  analysis  is  performed  for  a  one-dimensional  system,  we  discuss  in  the  conclusions  extensions 
of  this  analysis  to  multidimensional  systems. 


3.1  Properties  of  Green’s  Functions 

Consider  the  following  inhomogeneous  differential  equation: 

Cu  =  f  (19) 

where  both  the  solution  u  and  the  forcing  function  /  are  defined  on  the  interval  [0, 1]  with 
homogeneous  boundary  conditions.  Let  us  also  assume  that  the  spatial  differential  operator  C  is 
of  the  form 

£  =  E  (20) 

z£Z 

where  Z  is  a  subset  of  the  natural  numbers,  and  A^(r)  is  a  smooth  function  of  r,  i.e.  A2(r)  €  C°°. 

The  Green’s  function  integral  form  of  the  solution  can  be  written  as  the  following: 

u{r)=  f  G{r,s)f{s)ds  (21) 

Js 

where  G(r,  s)  is  the  Green’s  function  associated  with  C  and  is  itself  the  solution  to  the  following 
differential  equation. 

CGir,s)  =  5{r-s)  (22) 

where  C  is  applied  to  the  variable  r.  Note  that  we  use  r  to  denote  a  point  in  the  observation  or 
“output”  space  of  the  Green’s  function,  (e.g.  the  point  of  displacement  for  a  mechanical  system, 
while  we  use  s  to  denote  a  point  in  the  source  or  “input”  space  of  the  Green’s  function,  e.g.  the 
point  of  forcing). 

The  key  property  of  G  that  will  lead  us  to  its  compression  in  the  wavelet  basis  is  the  fact  that 
it  is  smooth  everywhere  except  along  its  diagonal  where  r  =  s,  i.e.  where  the  output  location 
coincides  with  the  input  location.  For  the  case  where  r  =  s,  from  eq.(22)  G  must  have  partial 
derivatives  of  some  order  that  are  singular.  One  way  to  see  this  is  to  observe  the  fact  that  the 
delta  function  on  the  right-hand  side  of  eq.(22)  implies  a  singularity  in  some  of  the  terms  on 
the  left-hand  side  of  eq.(22)  at  the  point  r  =  s.  From  eq.(20)  we  see  that  the  left-hand  side  of 
eq.(22)  comprises  a  sum,  weighted  by  smooth  functions,  of  partial  derivatives  of  G  with  respect 
to  components  of  r.  In  particular, 

^  A^(r)a;G(r,  s)  =  <^(r  -  s)  (23) 

z£Z 
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Let  r  +  1  denote  that  highest  derivative  in  the  summation  in  eq.(23).  We  then  rewrite  eq.(23), 
isolating  terms  involving  as  follows. 

\T+l{r)^QT^^G{r,s)  =  5{r  -  s)  +  e(r)  (24) 

where  e(r)  is  a  linear  combination,  weighted  by  smooth  functions,  of  derivatives  of  G  with  respect 
to  r  of  order  less  than  T+1.  Note  that  derivatives  of  G  with  respect  to  r  of  order  j  for  j  <  T  +  1 
must  be  strictly  smoother  than  the  T  +  1st  derivative  of  G.  Thus,  eq.(24)  shows  dJ'^^G  to 
be  smooth  for  all  r  except  at  r  =  s,  where  it  is  a  delta  function.  Note  that  this  makes  dj G 
discontinuous  at  r  =  s  and  G  G  C°°  for  r  ^  s.  We  state  this  fundamental  property  of  Green’s 
functions  as  follows. 


Property  3.1  Smoothness  Property  of  Green’s  Function  (I) 

Let  G(r;  s)  denote  the  Green’s  function  as  a  function  of  the  observation  point  r,  holding  the  source 
point  s  fixed.  Then,  G{r-,s)  G  C°°  for  r  ^  s  whereas  for  r  —  s,  G{r;s)  G  C^~^,  i.e.  djG{r;s)  is 
discontinuous  at  r  =  s. 


We  have  shown  that  by  fixing  s,  G  is  smooth  as  a  function  of  r  everywhere  except  at  r  =  s. 
However,  we  are  ultimately  interested  in  transforming  the  coordinates  of  the  Green’s  function 
operator  into  wavelet  basis  coordinates  in  both  input  and  output  spaces.  In  other  words,  in 
representing  our  integral  operator  we  need  to  integrate  G  with  respect  to  wavelet  basis  functions 
in  s  as  well  as  in  r.  Thus,  we  would  like  to  show  that  by  fixing  r,  G  is  smooth  as  a  function  of  s 
everywhere  except  at  s  =  r. 

In  order  to  describe  how  G  behaves  as  a  function  of  s,  holding  r  fixed,  we  invoke  the  reciprocity 
property  of  Green’s  functions  [10],  which  says  that  G{r,s)  =  G{s,r)  where  G(s,r)  is  the  Green’s 
function  associated  with  the  adjoint  of  £.  Whereas  before  we  analyzed  the  properties  of  £  to 
deduce  smoothness  properties  of  G  as  a  function  of  r,  we  now  analyze  properties  of  the  adjoint  of 
jC,  which  we  denote  as  £,  to  deduce  smoothness  properties  of  G  as  a  function  of  s.  In  particular, 
from  the  definition  of  £  in  eq.(20),  £  can  be  written  as  follows  [10]. 

£=S(-l)'a.'A.(s)  (25) 

z^Z 

The  Green’s  function  for  this  adjoint  operator,  G(s,r),  satisfies  the  following  equation. 

X:a;A.(3)G(»,r)=i(s-r)  (26) 

z^Z 


Given  our  assumption  of  smooth  weighting  functions  A2(r),  we  can  apply  the  same  argument 
to  eq.(26)  as  we  did  to  eq.(23).  It  then  follows  that  Property  3.1  holds  for  G  as  a  function  of 
s,  holding  r  fixed.  From  reciprocity,  we  can  then  conclude  that  Property  3.1  holds  for  G  as  a 
function  of  s,  holding  r  fixed.  In  particular, 
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Property  3.2  Smoothness  Property  of  Green’s  Function  (II)  Let  G(s;  r)  denote  the  Green’s  func¬ 
tion  as  a  function  of  the  source  point  s,  holding  the  observation  point  r  fixed.  Then,  G{s;  r)  €  C°° 
for  s-^r  whereas  for  s  =  r,  G(s;  r)  €  C^~^,  i.e.  dfG{s-,  r)  is  discontinuous  at  s  =  r. 

Finally,  we  note  that  for  self-adjoint  operators  the  reciprocity  property  of  Green’s  functions  states 
that  G{r,s)  =  G{s,r),  i.e.  the  Green’s  function  is  a  self-adjoint  kernel.  This  has  implications  on 
the  smoothness  properties  on  the  diagonal  of  G  in  that  drG{r,s)  =  dsG{s,r). 

Property  3.3  Smoothness  Property  of  Green’s  Function  (Self-Adjoint  Case)  Let  G{r,  s)  denote 
the  Green’s  function  of  a  self-adjoint  differential  operator.  Then,  G(r,s)  G  G°°  for  s  ^  r  whereas 
for  s  =  r,  dfd^G{r,s)  is  discontinuous,  whenever  a  +  ^  —  T  —  1. 

In  summary,  Green’s  functions  satisfying  eq.(19)  are  smooth  off  their  diagonal,  i.e.  for  r  ^  s, 
while  on  the  diagonal  they  are  (7^“^  where  T  -f-  1  is  the  highest  derivative  of  the  differential 
operator.  For  the  case  of  self-adjoint  differential  operators,  the  condition  on  the  diagonal  is  such 
that  mixed  partials  with  respect  to  r  and  s  are  discontinuous  when  the  sum  of  the  orders  of  the 
partials  is  T  —  1 . 


3.2  Structure  of  Wavelet  Transformed  Integral  Operator 


In  this  section  we  reveal  the  fundamental  structure  of  the  Green’s  function  integral  operator  in 
the  wavelet  basis.  We  begin  with  the  definition  of  the  operator  coefficients  in  the  wavelet  basis, 
assuming  a  complete,  orthonormal  wavelet  basis  defined  on  the  interval  [0, 1].  This  basis  set  is 
defined  as  follows. 

{$(2’"r-n)}  (27) 

for  all  integers  m,  n,  where  m  is  the  scale  index  and  n  is  the  translation  index. 

We  now  express  the  coefficients  of  the  integral  operator  in  the  basis  given  in  eq.(27)  as  follows. 

am,n,m,n  =  /  /  _  G{r,s)^{2‘"^r  -  n)$(2’”s  -  n)drds  (28) 

J2-^n  J2-^n 


Consider  the  following  equations  for  the  coefficients  of  the  transformed  functions  in  both  the 
input  and  output  spaces  of  the  operator. 


r2-^\ 
—  / 

,n  —  / 

J2-^n 


(n+l) 

^(2’"s  —  n)x{s)ds 

(29) 

(n+l) 

^(2^r  —  n)y{r)dr 

(30) 

Then,  the  following  equation  describes  the  integral  operator  mapping  between  these  functions  in 
the  wavelet  basis. 

'y  .  ^m.n  (^^) 


Vm^n  —  /  ^ 
7n,n 


no 
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Figure  2:  Matrix  representing  wavelet  transformed  integral  operator.  Dark  points  represent 
coefficients  corresponding  to  basis  functions  whose  support  include  the  diagonal  of  G.  A  8-tap 
Daubechies  filter  is  used  on  an  interval  with  256  uniformly  sampled  points  . 

In  performing  the  integrals  in  eq.(28)  there  are  two  cases  to  consider  based  on  the  bounds  of 
integration  with  respect  to  dr  and  the  bounds  of  integration  with  respect  to  ds,  which  we  denote 
as  follows. 

Clr  =  [2-’”n,2-’"(n-hl)]  (32) 

D,  =  [2-^n,2-^{n  +  l)]  (33) 

The  two  cases  to  consider  consist  of:  1)  flrDDs  =  0  and  2)  7^  0-  For  Case  1,  in  performing 

our  integration  we  can  consider  G{r,  s)  to  have  the  properties  for  the  case  of  r  /  s,  while  for 
Case  2,  in  performing  our  integration  we  must  consider  G{r,  s)  to  have  the  properties  for  the  case 
of  r  =  s.  An  interpretation  of  Case  1  is  that  we  are  considering  the  case  where  the  support  of 
the  wavelet  functions  for  the  observation  and  source  positions  do  not  overlap,  while  for  Case  2, 
we  are  considering  the  case  where  the  they  do  overlap. 

For  visual  purposes,  it  is  worthwhile  at  this  point  to  illustrate  the  structure  of  the  transformed 
operator  imposed  by  this  division  between  Case  1  and  Case  2.  In  particular.  Case  1  corresponds 
to  coefficients  whose  associated  basis  functions  include  the  diagonal  of  G  in  their  support.  The 
fact  that  basis  functions  that  include  the  diagonal  exist  at  multiple  scales  gives  rise  to  a  multiband 
structure,  where  each  band  corresponds  to  a  pair  of  basis  functions,  in  r  and  s,  whose  support 
differ  by  some  power  of  two.  The  width  of  each  band  is  simply  the  sum  of  the  supports  of  the 
basis  functions  in  r  and  s.  We  illustrate  this  multiple  band  structure  of  the  operator  in  Figure  2 
for  the  case  of  an  8-tap  Daubechies  filter  with  four  vanishing  moments. 

For  Case  1,  we  know  from  Properties  3.1  and  3.2  that  G(r,  s)  is  a  smooth  function  in  both  r 
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and  s.  Thus,  from  Taylor’s  theorem  we  can  expand  G(r,  s)  as  a  polynomial  in  r  and  s  up  to 
arbitrary  order.  Using  the  vanishing  moments  property  of  wavelets  in  eq.(9),  the  expression 
for  the  operator  coefficients  in  the  wavelet  basis  as  given  in  eq.(28)  can  be  reduced  to  a  double 
integral  of  the  remainder  term  of  the  Taylor  expansion  with  respect  to  the  wavelet  basis  functions 
in  r  and  s.  The  important  point  in  this  case  is  the  fact  that  for  wavelets  of  any  given  order 
vanishing  moments,  the  operator  coefficients  are  equal  to  a  double  integral  of  a  remainder  term 
corresponding  to  a  Taylor  expansion  of  that  order.  We  will  show  in  the  next  section  explicit  upper 
bounds  for  the  magnitude  of  these  coefficients.  In  particular,  this  will  show  that  for  Case  1  the 
magnitude  of  the  operator  coefficients  can  be  made  arbitrarily  small  given  wavelets  of  arbitrary 
order  vanishing  moments. 

For  Case  2,  we  know  from  Properties  3.1  and  3.2  that  for  r  =  s,  G(r,  s)  is  smooth  in  r  and  s  only 
to  a  certain  order  derivative  in  each  of  the  coordinates,  where  the  order  is  related  to  the  highest 
order  derivative  of  the  particular  differential  operator  C.  For  the  case  of  self-adjoint  £,  from 
Property  3.3  we  know  that  the  mixed  partials  with  respect  to  r  and  s  are  singular  when  the  sum 
of  the  orders  of  the  respective  partials  is  equal  to  the  highest  order  derivative  of  C.  Thus,  unlike 
in  Case  1,  we  cannot  use  a  Taylor  expansion  to  expand  G  up  to  arbitrary  order.  Furthermore, 
from  Corollary  2.1,  the  fact  that  G  is  smooth  only  up  to  some  order  in  each  coordinate  r  and  s 
implies  the  following.  For  each  coordinate  r  and  s,  if  n  —  1  is  the  order  of  the  highest  continuous 
partial  derivative  in  that  coordinate,  then  for  every  polynomial  of  order  i,  i  >  n,  G  in  that 
coordinate  cannot  be  approximated  by  that  polynomial  with  error  less  than 

From  Corollary  2.1,  we  can  make  the  stronger  claim  that  for  every  polynomial  of  order  i,  i  >  n, 
G  as  a  function  of  that  coordinate  cannot  be  approximated  by  a  polynomial  with  error  less  than 
i.e.  higher  order  approximations  of  G  have  an  error  that  is  independent  of  the  order  of 
the  approximation.  Thus,  in  contrast  to  Case  1,  where  G  could  be  expanded  to  arbitrarily  higher 
order  terms,  in  Case  2,  G  is  limited  to  polynomial  approximations  of  order  n  —  1  or  fewer  where  n 
reflects  the  structure  of  C.  Finally,  as  we  will  see  explicitly  in  the  next  section,  the  fact  that  we 
are  limited  in  the  order  of  the  polynomial  approximating  G  in  each  of  its  coordinates  limits  us 
in  the  number  of  vanishing  moments  that  can  be  used  to  bound  the  operator  coefficients  in  the 
wavelet  ba.sis.  This  is  a  key  distinction  in  discriminating  between  the  coefficients  in  Case  1  and 
the  coefficients  in  Case  2.  Roughly  speaking,  whereas  in  Case  1  the  magnitudes  of  the  coefficients 
can  be  made  ajbitrarily  small  with  higher  and  higher  order  vanishing  moments,  in  Case  2  they 
cannot. 


3.3  Bounds  on  Magnitude  of  Transformed  Operator  Coefficients 


We  now  derive  bounds  on  the  wavelet  transform  coefficients  for  a  Green’s  function  integral 
operator.  From  Properties  3.1  and  3.2,  we  know  G{r,s)  to  be  C°°  away  from  r  =  s  and  to  be 
T  —  1  times  differentiable  on  the  line  r  =  s^.  The  T’th  order  partial  derivatives  of  G(r,s)  are 

^Recall  that  differentiability  on  a  set  for  a  real  function  of  two  variables  requires  athat  the  partial  derivatives 
exist  and  are  continuous  on  the  set. 
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assumed  to  be  discontinuous  at  r  =  s  and  bounded  over  the  domain  of  interest. 

The  quadrupally  indexed  wavelet  coefficient,  a,  is  as  given  in  eq.(28),  where  ^  satisfies  the 
moment  property  stated  in  eq.(9),  i.e.  'J'  has  M  vanishing  moments.  Recall  that  the  bounds  of 
integration  with  respect  to  dr  and  ds,  respectively,  in  eq.(28)  are  given  by 

Ur  =  [2-™n,2-'”(n  +  l)],  and 
a  =  [2-’^n,2-^(h  +  l)] 

which  reflect  the  support  of  the  wavelet  basis  function  in  the  r  and  s  directions,  respectively. 
We  derive  bounds  for  the  two  cases  identified  in  the  previous  section:  Case  1)  fir  H  f)*  =  0  and 
Case  2)  fir  H  fis  ^  0.  The  bounds  rely  upon  the  smoothness  properties  of  G{r,  s)  and  use  the 
multidimensional  Taylor  series  expansion  along  with  the  theorem  and  corollaries  introduced  in 
Section  2.2. 


3.3.1  Case  1:  Non-intersecting  wavelet  support 

When  the  support  of  the  wavelet  bases  in  r  and  s  not  overlap,  the  Green’s  function  in  the  integral 
in  eq.(28)  will  be  C°°  over  the  entire  domain  of  the  integral.  We  derive  an  upper  bound  for  the 
magnitude  of  the  wavelet  coefficient  as  a  function  of  scale  as  follows. 

Since  G{r,s)  is  C°°  for  r  and  s  of  interest,  we  can  expand  it  in  a  two-dimensional  Taylor  series 
about  a  point,  (ro,  5o))  to  arbitrary  order.  For  convenience,  choose  the  expansion  point,  (ro,  sq) 
to  be  (2~’"n,  2~’"n),  the  endpoint  of  the  wavelet  support.  We  choose  to  expand  the  first  2M  —  1 
polynomial  terms  in  the  Taylor  series  so  that  we  can  apply  the  moment  property  of  the  wavelet 
to  annihilate  the  low  order  behavior  of  the  Green’s  function.  The  Taylor  series  with  remainder 
is  given  by 

2M-1  j  /  Q  q\3 

G{r,s)  =  G{ro,so)+  Xj  G'('’’^)l'-=ro,s=*o  +  ^2M(r,s)  (34) 

The  above  expansion  involves  no  approximation  when  the  remainder  is  included. 

The  remainder  is  given  by 

R2M+i(r,s)  =  Jo  (('■  “  +  (*  -  G(ro+t{r-ro),so+t(s-so))dt 

(35) 

By  substituting  the  above  expansion  for  G  into  eq.(28),  applying  the  moment  property  of  the 
wavelet  basis  annihilates  all  of  the  polynomial  terms  in  the  expansion,  leaving  for  a,  without 
approximation,  the  following. 

oim,n,m,n  =  -  n)$(2”’s  -  n)i?2M(r,  s)drds  (36) 
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We  can  bound  a  by  applying  the  Cauchy-Schwartz  inequality  to  the  above  equation  giving 
<  [  I  '$^{2'^r  —  n)'^^{2'^s  —  n)drds  /  /  Rlj^{r,  s)drds  (37) 

J2-”'n  J2-'^n  J2-”'n  J2~*^n 

The  orthonormal  property  of  the  wavelet  basis,  implies  that  the  first  double  integral  is  unity. 
To  construct  the  bound,  we  note  that 


a  a  \  2M  a2M 

i=0 


'5r  '  dr'ds^M-i 

where  (3i  axe  positive  integer  coefficients  of  a  binomial  expansion.  Therefore,  we  can  write 

^2m(’'5^)  =  (2M  -  1)!2  ~  ■“  ’0.J 


(38) 


e, 


Jo 


_1  52^ 


dr'‘ds^^ 


—■G{rQ  +  t{r  -  To),  5o  +  <(5  -  SQ))dt 


(39) 


Furthermore,  we  can  use  the  inequality^ 


/2M  \  2  2M 

^a,  <(2M  +  l)j:a? 

Vt=0  /  i=0 


to  show  that 


(40) 


We  bound  the  integral  in  parenthesis  as  follows: 

=  (£(1  -  -  ^0))^^) 

/•I  rl  (  \2 

<  J  {\  —  tY^~^dtj  (  ^--T^-^^j^G(ro  +  t(r  —  ro),5o  +  t(5  —  5o))dt  j  dt  (Cauchy-Schwartz) 


< 


B? 


(4M  -  1) 
where  the  bound  Bf  is  defined  to  be 


(41) 


/  a2M  \  2 

=  sup  s)dtj  r  ^  s 


^The  Cauchy-Schwartz  Inequality  for  sequences  is  ^  5^i=i  X^i=i  *  Setting  6,-  =  1  yields 
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Returning  to  eq.(40),  we  can  write 


2M 


AM-2i 


(42) 


Finally,  returning  to  eq.(37)  and  substituting  eq.(42)  gives  for  a  bound  on  a 

’(n+l)  />2-’»(n+l)  2Af 


a2< 


(2M  +  1) 


(4M  -  1)(2M  -  1)!2 


f  f  ^  -  ^o)  ‘(^  -  ^o) 

J2-’^n  J2-*^n  ^ 


4M-2t 


drds 


Evaluating  the  integrals  in  r  and  s  gives  for  the  bound 

2  (2M  +  1)  ^ _ pjBf _ 

-  (4M  -  1)(2M  -  1)!^  ^  (2i  +  1)(4M  -  2t  +  1) 


(43) 


(44) 


Defining  a  single  coefficient, 

^ _ (2M  +  mBf _ 

(4M  -  1)(2M  -  l)!2(2z  +  1)(4M  -  2i  +  1) 

highlights  the  scale  dependence  of  the  bound 

2M 

^m,n, m,n  ^  ^ 

t=0 


(45) 


(46) 


The  above  equation  bounds  the  wavelet  coefficient  a  function  of  scale,  m  and  m,  for  Case  1. 
The  only  parameter  determining  the  rate  of  decay  with  scale  is  the  wavelet  moment  order,  M. 
Note  that  the  bound  can  be  arbitrarily  decreased  by  increasing  the  number  of  vanishing  moments, 
M.  Numerical  confirmation  of  this  bound  will  be  shown  in  Section  4. 


3.3.2  Case  2:  Intersecting  wavelet  support 

We  now  turn  to  the  case  where  the  intersection  of  Dr  and  D*  is  not  empty.  Here,  the  integral 
expression  in  eq.(28)  contains  the  point  r  =  s  for  which  the  Green’s  function  is  discontinuous 
in  the  T’th  order  partial  derivatives.  We  again  derive  a  scale  dependent  bound  for  the  wavelet 
coefficient  in  this  case  that  closely  follows  the  derivation  in  Case  1.  We  show  below  that  the 
bound  in  this  case  will  depend  upon  the  relationship  between  the  “depth”  of  the  singularity  in 
G{r,s),  T,  and  the  moment  property  of  the  wavelet  basis,  M.  When  2M  <  T  —  1,  the  wavelet 
basis  does  not  annihilate  enough  of  the  low  order  behavior  of  the  Green’s  function  to  uncover 
the  singularity,  and  the  bound  is  the  same  as  that  derived  in  Case  1.  However,  if  2M  >  T  —  1, 
the  wavelet  strips  away  enough  of  the  low  order  behavior  of  G{r,  s)  to  expose  its  singularity,  and 
the  bound  on  the  wavelet  coefficient  differs  from  Case  1.  The  derivation  for  Case  2  proceeds  as 
follows. 
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As  in  Case  1,  we  expand  G(r,  s)  about  a  point  (ro,so)  =  (2'"n,  2'"n).  Given  that  G{r,s)  is 
at  r  =  s,  we  can  expand  G  in  the  first  T  —  1  terms  of  its  Taylor  series 

1  /  d  d  \  ^ 

G(r,s)  =  G(ro,so)  +  Z)  +  (^  -  G{r,s)\r=ro,s=so  +  RT{r,s)  (47) 

with  the  remainder  given  by 

Rrir, s)  =  ,y  _  j,|  /  j  +  *('■“  +  *(«  -  *o))<ft 

(48) 

Since  we  assume  that  G{r,s)  has  bounded,  discontinuous  derivatives  at  r  =  s,  the  remainder 
above  is  well  defined.  Note  that  the  expansion  we  are  performing  in  eq.(47)  is  essentially  a  one¬ 
dimensional  Taylor  expansion  of  G  along  the  line  connecting  (ro,so)  and  (r,s);  i.e.  (ro  -t-  t{r  - 
f’o),so  +  t{s  —  So)).  The  expansion  is  performed  with  respect  to  the  variable  t,  which  represents 
the  extent  along  this  line.  Thus,  by  Theorem  2.1  we  know  that  if  we  wish  to  approximate  G  by  a 
polynomial  of  order  T  —  1  along  this  line  with  0{t^)  error,  the  Taylor  expansion  in  eq.(47)  does 
precisely  this.  Furthermore,  by  Corollary  2.2  of  Section  2.2,  if  we  wish  to  approximate  G  along 
(ro  4-  tr.  So  +  ts)  by  a  polynomial  of  order  T  or  higher  the  error  is  greater  than  Thus, 

given  that  the  Green’s  function  is  discontinuous  in  the  T’th  order  partial  derivatives  for  r  =  s, 
eq.(47)  represents  the  highest  order  expansion  of  G  with  0{t^)  error;  higher  order  expansions  do 
not  result  in  lower  order  error. 

Now,  when  we  substitute  eq.(48)  into  the  integral  equation  for  the  wavelet  coefficient,  eq.(28), 
we  make  the  following  observation.  If  the  moment  order  M  is  such  that  2M  <  T  —  1,  then  the 
wavelet  bases  will  not  annihilate  all  of  the  polynomial  terms  in  the  expansion  for  G,  eq.(47).  The 
remaining  polynomial  terms  will  dominate  the  wavelet  decay  as  a  function  of  scale,  precisely  as  in 
Case  1,  and  the  bound  on  the  wavelet  coefficients  as  a  function  of  scale  will  be  exactly  that  given 
in  eq.(46).  Conversely,  if  2M  >  T  —  1,  the  wavelet  bases  will  annihilate  all  of  the  polynomial 
terms  in  the  expansion  for  G,  and  the  behavior  of  the  wavelet  coefficient  will  be  determined  by 
RT(r,s).  We  turn  now  to  deriving  the  bound  in  this  case,  2M  >  T  —  1. 

The  derivation  of  the  bound  follows  that  of  Case  1.  Substituting  the  expansion  of  G  in  eq.(48) 
into  eq.(28)  and  using  the  moment  property  of  the  wavelet  bases  with  2M  >  T  —  1  leaves  for  a 

o;m,n,m,n  =  /  /  ’I'(2’”r  -  n)'^ (2™s  -  n)/?7-(r,  s)drds  (49) 


Using  the  Cauchy-Schwartz  inequality  and  the  orthonormal  property  of  the  wavelet  bases  gives 


As  before,  we  can  bound  as  follows 

^  =  (r  -  1)!2  ~  ~ 

'T 

0i  =  Jo^^~  0  +  t{r  -  To),  So  +  tis  -  So))dt 

^  -  {T-iy^  E  ~  (0.f 

^  -  {2T  -^Kr  -  1)!^  ~  ~ 

where  /3,  are  positive  integer  coefficients  of  a  binomial  expansion  and  Bf  is  given  by 

Bf  =  Bup 


(51) 


(52) 

(53) 


(54) 


Returning  to  eq.(50)  and  substituting  eq.(54)  bounds  a  as 

r2-"*(n+l) ,  T 


< 


(T+l) 


L  L  E^^mr-rons-so) 

J2-’"n  J2-*^n  ^ 


(2T-l)(r-l)!2 


2T-2i 


drds 


(55) 


Finally,  evaluating  the  integrals  in  r  and  s  gives 

2  ^  (r  +  1)  ^ _ PiBj _ «-(2i+l)m«-(2T-2j+2)m 

-  {2T-  1){T  -  1)!2  {2i  +  l)(2r  -  2i  +  1) 

As  before,  we  can  highlight  the  scale  dependence  of  the  bound  by  defining 


-  ^ _ {T  +  l)Mf _ 

(2r-l)(r-l)!2(2i  +  l)(2T-2i  +  l) 


to  give  for  the  bound 


T 

^2  ^  ^ -^.2“(2«+l)m2~(2T-2i+2)m 

t=0 


(56) 


(57) 

(58) 


The  above  bound  gives  the  decay  of  the  wavelet  coefficient  with  scale  when  the  support  of  the 
wavelet  includes  the  singularity  in  the  Green’s  function.  It  depends  upon  the  smoothness  of  the 
Green’s  function,  T,  and  the  scales  of  the  wavelet,  m  and  m  and  is  applicable  when  2M  >  T  —  1. 
Note  that  the  form  of  this  bound  is  identical  to  that  of  CaseT  with  M  in  eq.(46)  being  replaced 
by  T.  This  highlights  the  fact  that  whereas  for  Case  1  the  bound  on  the  coefficients  can  be  made 
arbitrarily  small  by  increasing  the  number  of  vanishing  moments  of  the  wavelet,  for  Case  2  the 
bound  for  2M  >  T  —  1  is  independent  of  M.  Numerical  confirmation  of  this  bound  will  be  shown 
in  Section  4. 
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3.4  Compression  Ratio 


Finally,  we  conclude  our  analysis  of  the  wavelet  representation  of  the  Green’s  function  integral 
operator  with  a  discussion  of  the  degree  of  compression  one  achieves  by  approximating  the 
integral  operator  by  its  significant  wavelet  coefficients.  The  bounds  in  eq.’s(46,58)  show  that  by 
using  wavelets  with  an  arbitrary  number  of  vanishing  moments,  we  can  make  the  coefficients  for 
Case  1  drop  off  substantially  faster  away  from  the  diagonal  than  coefficients  for  Case  2.  This  leads 
naturally  to  the  idea  of  approximating  the  operator  by  zeroing  out  coefficients  corresponding  to 
Case  1,  while  retaining  coefficients  corresponding  to  Case  2.  Such  an  approximation  gives  rise 
to  an  operator  whose  multiband  structure  is  as  illustrated  in  Figure  2.  We  now  derive  a  formula 
for  the  number  of  coefficients  remaining  in  the  approximation  as  a  function  of  the  number  of 
samples  in  the  interval,  2M  the  filter  length,  and  L,  the  number  of  scales. 

Consider  N  samples  on  the  interval  [0, 1],  giving  rise  to  an  TV  x  matrix  representation  of  the 
Green’s  function  integral  operator.  We  let  K  denote  the  number  of  significant  coefficients  in  the 
wavelet  representation  of  the  integral  operator,  i.e.  the  total  number  of  coefficients  at  all  scales 
lying  within  the  bands  representing  the  case  where  0.  We  denote  the  number  of  scales 

as  Ir  =  logj  —  1.  Recall  that  the  filter  length  is  2M.  Let  W  denote  the  width  of  the  band  at 
scales  m,  fn.  Based  on  the  fact  that  W  represents  the  coefficients  at  particular  scales  m  and  fn, 
where  fir  H  fls  0,  we  get  the  following  formula  for  W. 

W  =  M  +  M2l’”-^l  -  1  (59) 


Summing  all  the  bands  in  the  transformed  operator  matrix,  we  get 

L  m  L 

K  =  ^2^  -  E  ^2”* 

m=0  m=0  m=0 

L  m  L 

=  2  ^  ^((M-l)2^+M2’”)-  52(2M-1)2^ 

771=0  7n=0  771=0 

Note  that 

L  m  _ 

^  ^  2”^  =  2(2^+^  -  1)  -  (L  +  1) 

771=0  771=0 

We  also  need  to  compute  the  following  sum. 

L  m  L 

E  E2”=  E2”(m  +  1) 

771=0  771=0  771=0 


The  sum,  can  be  bounded  above  by 

/  2Hdt  =  7(2^+*  (L  -  7  +  1)  +  7) 

Jo 


(60) 


(61) 


(62) 


(63) 
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where  7  =  7^- 
Thus, 

K  <  2{(M-1)(2(2^+'-1)-(L  +  1))  +  M(7(2^+HL-7  +  1)  +  7)  +  (2^+1-1))} 

-  (2M  -  1)(2^+^  -  1) 

<  (6M  -  A){N  -  1)  +  7(T  -  7  +  l)2MiV  -  (2M  -  2){L  +  1)  +  2M72 

-  (2M  -  l)N  (64) 


For  large  N, 


K  =  0({2M(7(I  -  7  +  1)  +  2)  -  3}iV) 

=  C>((4M  +  27ZM  +  27(1 -7)M-3)iV) 


(65) 


From  eq.(65)  we  see  that  the  number  of  coefficients  in  the  approximation  to  the  operator  is  of 
order  linear  in  as  opposed  to  as  in  the  case  of  the  full  operator.  The  number  of  significant 
coefficients  is  also  dependent  on  filter  length  and  the  number  of  scales.  Clearly,  this  approxima¬ 
tion  becomes  advantageous  rather  quickly  for  large  N.  Moreover,  by  using  wavelets  with  higher 
order  moments  one  can  arbitrarily  increase  the  magnitude  of  the  ratio  of  the  coefficients  on  the 
bands  to  the  coefficients  off  the  bands,  thus  decreasing  the  approximation  error.  This  is  done, 
however,  at  the  cost  of  increasing  M  and  so  there  is  a  natural  tradeoff  between  approximation 
error  and  number  of  coefficients. 


4  Numerical  Examples 


The  results  of  the  prior  section  allow  us  to  bound  the  wavelet  coefficient  of  the  Green’s  function 
subject  to  the  two-dimensional  wavelet  transform  defined  in  eq.(  28).  When  the  support  of  the 
wavelets  in  r  and  s  do  not  overlap,  the  bound  on  the  wavelet  coefficient  is  given  by  eq.(  46). 
This  bounds  wavelet  coefficients  off  of  diagonal  and  subdiagonal  blocks  in  the  wavelet  transform 
space. 

When  the  support  of  the  r  and  s  wavelets  overlap  and  2M  <  T  —  1,  so  that  the  wavelet 
moment  property  is  insufficient  to  annihilate  enough  low  order  behavior  of  G(r,  s)  to  uncover 
the  singularity,  eq.(  46)  still  bounds  the  wavelet  coefficients.  However,  when  2M  >  T  —  1,  the 
wavelet  strips  away  the  low  order  behavior  of  the  Green’s  function,  up  to  and  including  the 
highest  order  continuous  derivative,  and  is  bounded  by  eq.(  58).  A  numerical  example  confirms 
these  conclusions. 

For  our  example,  we  use  an  Euler-Bernoulli  beam,  whose  equation  of  motion  is  given  by 

Wrrrr  +  =  -p{r)  T  €  [0,  1] 
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Figure  3:  A  plot  of  the  magnitude  of  the  Green’s  function  for  our  Euler-Bernoulli  beam  with 
clamped  ends  as  a  function  of  the  field  point,  r  and  the  source  point,  s,  and  dimensionless 
frequency,  u)  =  215.  The  scale  is  given  by  the  colorbar  on  the  right. 

where  w  is  the  dimensionless  vertical  displacement  of  the  beam,  uj  is  dimensionless  frequency 
and  p  represents  the  forcing  along  the  beam.  We  use  clamped  boundary  conditions 

10  =  Wr  =  0  for  r  =  0, 1. 


The  Green’s  function  in  this  case  satisfies 

G(r,  s)rrrr  +  s)  =  -5{r  -  s), 

along  with  the  boundary  conditions  above.  Solution  for  the  Green’s  function  is  possible  in 
closed  form.  Figure  3  shows  the  magnitude  of  the  Green’s  function  at  a  dimensionless  frequency, 
(jj  =  215,  used  for  this  example. 

The  Green’s  function  for  the  Euler-Bernoulli  beam  above  has  T  =  3,  so  that  G(r,  s)  is  (7°° 
everywhere  except  on  r  =  s  where  it  is  C^.  Off  of  the  bajids  in  the  wavelet  representation  of  the 
operator,  i.e.  Case  1,  we  expect  the  wavelet  coefficients  to  decay  in  scale  as 

2M 

C?  -  -  <  ^..y^.2“(2i+l)m2-(4Af-2»+l)m^ 

1=0 

where  M  is  the  number  of  vanishing  moments  of  the  wavelet.  In  our  examples  we  use  the 
Daubechies  wavelets. 
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Figure  4:  Decay  of  numerically  calculated  wavelet  coefficient  a,  with  scale  along  diagonal  blocks 
of  the  wavelet  transform,  m  =  fh.  Slope  is  predicted  to  go  as  — (4M  +  2).  The  kink  in  the  M  =  5 
curve  is  due  to  the  fact  that  the  m  =  2  point  is  smaller  than  machine  precision. 


We  test  our  bound  by  computing  a  numerically  and  finding  the  supremum  of  a  off  of  the  diagonals 
and  subdiagonals  for  each  scale,  m  and  fh.  We  can  plot  our  results  compactly  by  noting  that 
along  any  set  of  scales,  m  =  m  +  fc,  A:  an  integer,  the  bound  looks  like 

log2(Q:^)  <  —  (4M  +  2)m  +  C. 


Therefore,  we  expect  21og2(Q:)  to  be  linear  with  slope  — (4M+2)  as  a  function  of  m  along  diagonal 
blocks.  Figure  4  shows  the  decay  of  the  supremum  of  log2(a)  along  the  blocks  m  =  fh  when  the 
wavelet  bases  do  not  intersect.  The  lines’  intersections  with  the  m  =  0  axis  have  been  normalized 
to  zero  to  show  their  slopes.  Shown  are  the  decay  rates  for  various  wavelet  bases  having  different 
numbers  of  vanishing  moments.  As  more  moments  are  annihilated  by  the  wavelet,  M  increasing, 
we  expect  faster  decay  of  the  wavelet  coefficient,  with  slopes  given  by  —  (4M  +  2).  Figure  4 
shows  that  this  is  indeed  the  case.  The  slopes  in  the  figure  show  good  agreement  with  the 
theoretical  bound.  We  believe  that  departures  from  the  theoretical  values  are  due  to  the  discrete 
implementation  of  the  wavelet  transform.  Note  that  the  last  point  in  the  M  =  5  curve  is  smaller 
than  machine  precision  and  is,  thus,  unreliable.  In  Figure  5,  we  show  the  case  where  m  =  fh+1. 
Again,  the  decay  as  a  function  of  scale  is  well  predicted  by  our  bound. 

Next,  we  test  the  bound  for  the  wavelet  coefficients  when  they  are  on  the  diagonals  of  the  wavelet 
transform  blocks.  Recall  that  we  expect  the  Case  1  bound  if  2M  <  T  —  1  and  the  Case  Two 
bound  if  2M  >  T  —  1.  For  the  Euler-Bernoulli  beam,  T  =  3,  so  that  only  the  M  =  1  wavelet 
(the  Haar  wavelet),  will  give  the  Case  1  bound.  Conversely,  when  2M  >  T  —  1,  we  derived  a 
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Figure  5:  Decay  of  numerically  calculated  wavelet  coefficient,  a,  as  a  function  of  scale  for, 
fh  =  m  +  \.  Theoretically  predicted  slopes  are  given  by  —{AM  +  2). 

bound  whose  scale  dependence  is  given  by 

T 

q,2  <;  ^ -^2“(2«+l)m2-(2T-2i+2)m^ 

«=0 

Again,  consider  the  case  where  m  =  fh  +  k,  with  k  an  integer.  For  M  =  1,  the  Haar  wavelet,  we 
expect  the  Case  2  bound  to  hold  on  the  main  diagonal  so  that 

log2(Q:^)  =  -6m  +  Co, 

along  the  diagonal.  However,  for  M  >  1,  we  expect  the  Case  Two  bound  to  hold  so  that 

log2(a^)  =  — (2T  +  2)m  =  — 8  +  Ci. 

Figure  6  shows  the  decay  of  numerically  obtained  wavelet  coefficients,  supremum  of  a,  along 
the  main  diagonal  m  =  m,  as  a  function  of  scale  for  wavelets  with  M  —  1,2,3  and  4.  As 
predicted,  the  slope  of  the  decay  for  M  =  1  is  —6  while  for  all  M  >  1  the  slope  is  —8.  Likewise, 
Figure  7  shows  the  decay  of  the  wavelet  coefficient  as  a  function  of  m  along  the  first  subdiagonal, 
fh  —  m  +  1-  As  before,  we  expect  the  slope  of  the  line  to  be  —6  for  M  =  1  and  —8  for  M  >  1. 
Numerical  agreement  with  the  predicted  bound  is  quite  good,  as  shown  in  the  figure. 

Finally,  we  show  in  Figure  8  the  wavelet  transformed  operator  using  an  8-tap  Daubechies  wavelet. 
Note  the  multiband  structure  reflecting  the  ability  of  the  wavelet  to  discern  the  singular  behavior 
of  the  Green’s  function  along  its  diagonal.  Figure  9  shows  the  approximated  operator  derived  by 
keeping  only  coefficients  whose  magnitudes  are  greater  than  le-8  times  the  norm  of  the  operator. 
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Figure  6:  Decay  of  numerically  calculated  wavelet  coefficient,  a,  as  a  function  of  scale  along  the 
main  diagonal.  Our  theory  predicts  a  slope  of  -6  for  M  =  1  and  a  slope  of  —8  for  M  >  0.  The 
above  plot  shows  excellent  agreement  with  the  theoretical  bounds. 


Figure  7:  Decay  of  numerically  calculated  wavelet  coefficient,  a,  as  a  function  of  scale  along  the 
first  sub-diagonal,  m  =  m  -1- 1.  Our  theory  predicts  a  slope  of  -6  for  M  =  1  and  a  slope  of  —8 
for  M  >  0.  The  above  plot  shows  excellent  agreement  with  the  theoretical  bounds. 
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Index  for  m  and  n 

Figure  8:  Wavelet  transformed  operator  using  8-tap  Daubechies  wavelet.  A  dB  intensity  scale  is 
used. 


Index  for  m  and  n 


Figure  9;  Thresholded  values  (shown  in  black)  of  wavelet  transformed  operator  using  8-tap 
Daubechies  wavelet.  Values  greater  than  le-8  times  the  norm  of  the  operator  are  shown. 
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5  Conclusions 


In  this  paper  we  have  analyzed  the  use  of  wavelets  with  vanishing  moments  in  representing 
Green’s  function  integral  operators  associated  with  linear  differential  equations.  We  have  shown 
that  this  representation  is  sparse,  manifesting  itself  in  a  multiband  structure,  and  that  this 
structure  follows  from  the  moment  property  of  the  wavelet  and  the  structure  of  the  singularity  of 
the  Green’s  function  along  its  diagonal.  In  particular,  we  show  that  Green’s  functions  are  smooth 
off  its  diagonal  and  limited  in  its  smoothness  on  its  diagonal  in  the  sense  that  at  some  order  its 
partial  derivatives  are  discontinuous,  where  this  order  corresponds  to  the  highest  derivative  of  the 
associated  differential  operator.  We  then  show  that  in  the  representation  of  the  integral  operator 
in  the  wavelet  basis,  coefficients  corresponding  to  wavelet  basis  functions  whose  support  include 
the  diagonal  of  the  Green’s  function  behave  much  differently  from  coefficients  corresponding  to 
basis  functions  whose  support  does  not  include  the  diagonal.  Namely,  coefficients  that  do  not 
include  diagonal  information  can  be  made  arbitrarily  small  as  the  number  of  vanishing  moments  of 
the  wavelet  increase,  whereas  coefficients  that  do  include  diagonal  information  cannot.  Moreover, 
the  bound  for  the  magnitude  of  these  latter  coefficients  is  a  function  of  the  highest  derivative  of 
the  associated  differential  operator. 

The  results  of  this  analysis  show  the  wavelet  representation  of  the  integral  operator  to  have 
multiple  bands,  each  of  which  represents  coefficients  corresponding  to  wavelet  basis  functions  of 
a  certain  scale  whose  support  includes  the  diagonal  of  the  Green’s  function.  Our  bounds  support 
the  idea  of  approximating  the  operator  by  zeroing  out  coefficients  off  of  the  bands,  since  these 
coefficients  drop  off  substantially  faster  away  from  the  diagonal  than  coefficients  on  the  bands. 
The  compression  ratio  one  gets  from  doing  this  is  reflected  by  the  fact  that  the  approximated 
operator  has  0((4M  +  2'jLM  +  27(1  —  7)M  —  2)N)  coefficients  as  opposed  to  N'^.  Finally, 
in  Section  4  numerical  results  show  good  agreement  with  predicted  results  for  the  case  of  an 
Euler-Bernoulli  beam  whose  behavior  is  governed  by  a  fourth  order  differential  equation. 

We  are  currently  extending  the  results  of  this  paper  to  the  case  of  PDF’s.  In  this  case  the 
smoothness  characterization  of  Green’s  functions  is  more  subtle  than  in  the  case  of  one  dimen¬ 
sional  problems.  In  particular,  the  relationship  between  the  derivatives  of  the  partial  differential 
operator  and  the  smoothness  of  its  associated  Green’s  function  along  its  diagonal  needs  to  be 
more  fully  understood.  We  believe,  however,  that  results  analogous  to  those  in  this  paper  can 
be  derived  for  the  PDF  case. 

Acknowledgements 

This  work  is  being  performed  under  an  ARPA  program  sponsored  by  the  Office  of  Naval  Research 
under  Contract  No.  N00014-94-C-0176. 


125 


References 


[1]  G.  Beylkin,  R.  Coifman,  and  V.  Rokhlin,  “Fast  Wavelet  Transforms  and  Numerical  Algo¬ 
rithms  I,”  Communications  on  Pure  and  Applied  Mathematics,  Vol.  XLIV,  pp.  141-183,  1991. 

[2]  K.  Chou,  G.  Guthart,  and  D.  Flamm,  “A  Multiscale  Approach  to  the  Control  of  Smart 
Materials,”  SPIE  Conference  on  Smart  Structures  and  Materials,  San  Diego,  March  1995. 

[3]  W.  Dahmen,  S.  Prossdorf,  and  R.  Schneider,  “Multiscale  Methods  for  Pseudo-differential 
Equations  on  Smooth  Closed  Manifolds,”  in  Wavelets:  Theory,  Algorithms,  and  Applications, 
C.  Chui,  L.  Montefusco,  L.  Puccio  (editors).  Academic  Press,  1994,  pp.385-424. 

[4]  I.  Daubechies, “Orthonormal  bases  of  compactly  supported  wavelets”.  Communications  on 
Pure  and  Applied  Mathematics,  Vol.  XLI,  pp.  909-996,  1988. 

[5]  S.  Jaifard, “Wavelet  Methods  for  Fast  Resolution  of  Elliptic  Problems,”  SIAM  Journal  on 
Numerical  Analysis,  Vol.  29,  No.  4,  pp.  pp.  965-986,  August  1992. 

[6]  S.G.  Mallat,  “A  Theory  for  Multiresolution  Signal  Decomposition:  The  Wavelet  Represen¬ 
tation”,  IEEE  Transactions  on  Pattern  Anal,  and  Mach.  Intel.,  Vol.  PAMI-11,  July  1989,  pp. 
674-693. 

[7]  S.  Mallat  and  W.L.  Hwang,  “Singularity  Detection  and  Processing  with  Wavelets,”  IEEE 
Transactions  on  Information  Theory:  Special  Issue  on  Wavelet  Transforms  and  Multiresolution 
Signal  Analysis,  Vol.  38,  No.  2,  March  1992,  pp.  617-643.. 

[8]  C.  Meneveau,  “Analysis  of  Turbulence  in  the  Orthonormal  Wavelet  Representation,”  J.  Fluid 
Mechanics,  Vol.  232,  1991,  pp.  469-520. 

[9]  Y.  Meyer,  Ondelettes  et  Operateurs,  Paris,  France:  Herman,  1990. 

[10]  P.  Morse  and  H.  Feshbach,  Methods  of  Theoretical  Physics:  Part  I,  McGraw-Hill,  1953. 

A  Appendix 


Lemma  A.l  For  n  a  positive  integer,  if  a  function  f(x)  is  Lipschitz  n  1  at  xq,  then  it  is 
Lipschitz  n  at  xq  and  Pn{h)  =  P'n-i{h)  +  phA ,  where  Pn{h)  and  P„_i(/i)  are  the  polynomials  in 
the  Lipschitz  n  +  1  and  n,  respectively,  definitions  for  f  at  the  point  Xq. 

Proof 

We  first  assume  f{x)  is  Lipschitz  n-(-l  at  xq.  From  Definition  2.1  this  means  there  exist  constants 
A  and  ho  >  0,  and  a  polynomial  of  order  n,  Pn{x),  such  that  for  0  <  <  ho 

l/(xo  +  h)-Pn(h)|<  (66) 
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Let  Pn{h)  =  Pn-\{h)  +  ph'^,  where  Pn-i{h)  represents  the  first  n  terms  of  the  polynomial  P„(h). 
Substituting  this  into  eq.(66)  and  using  the  triangle  inequality  we  get 

|/(xo  +  h)-  Pn-i{h)  \  -  |p|/i"  <  (67) 

which  can  be  rewritten  as 

|/(xo  +  h)-  Pn-i{h)\  <  (IpI  +  Ah)h'"  (68) 

If  we  let  A  =  \p\-\-  Aho,  then 

|/(xo  +  h)-P„_i(/i)|  <>lV  (69) 

By  definition  /  is  Lipschitz  n  at  xq,  where  the  polynomial  in  the  definition  is  in  fact  P„_i(h), 
i.e.  the  first  n  terms  of  the  polynomial  in  the  definition  /’s  being  Lipschitz  n  +  1  at  xq. 

□ 


Theorem  A.l  For  n  a  positive  integer,  if  /(x)  is  Lipschitz  n  +  1  at  xq  and  /(x)  is  (7"  for 
0  <  |x  —  xo|  <  7,  7  >  0;  then  /(x)  is  n  times  differentiable  at  xq  and  the  polynomial  Pn{h)  is 
equal  to  the  first  n  +  1  terms  of  the  Taylor  series  expansion  of  f[x)  about  xq. 


Proof 


We  use  induction  on  n  to  prove  this  result.  For  n  =  1,  /(x)  is  Lipschitz  2  at  xq.  By  letting  h  =  0 
in  eq.(12),  we  get  the  following. 

|/(xo)  -  Pi(0)|  <  0  (70) 

Thus,  Pi(0)  =  /(xo).  By  using  this  fact  and  Definition  2.1,  we  know  that  there  exist  constants 
A  and  ho  >  0,  and  a  polynomial  of  order  1,  Pi(x),  such  that  for  0  <  h  <  /iq 

\f{xo  + h)  -  f{xQ)  -  pih\<  Ah'^  (71) 


where  pi  is  the  coefficient  of  the  first  order  term  of  Pi{h).  Dividing  both  sides  of  eq.(71)  by  h, 
we  get 


For  every  c  >  0,  we  can  choose  h  such  that  the  left-hand  side  of  eq.(71)  is  bounded  by  e  whenever 
h  <h.  Thus,  the  first  derivative  of  /  at  xq  exists  and  is  in  fact  equal  to  pi. 


/(xo  +  h)-  /(xo) 


—  Pi  I  <  Ah 


(72) 


Our  induction  hypothesis  is  stated  as  follows.  For  n  —  1  a  positive  integer,  if  /(x)  is  Lipschitz  n 
at  Xo  and  /(x)  is  (7"“^  for  0  <  |x  -  xo|  <  7,  7  >  0,  then  /(x)  is  n  —  1  times  differentiable  at  xq 
and  the  polynomial  Pn-i{h)  is  equal  to  the  first  n  terms  of  the  Taylor  series  expansion  of  /(x) 
about  Xq. 
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We  now  show  our  induction  hypothesis  is  true  for  n.  We  assume  in  our  theorem  statement  that 
/(x)  is  Lipschitz  n  +  1  at  xq.  From  Lemma  A.l,  we  know  that  /(x)  is  Lipschitz  n  at  Xq,  and 
thus,  from  our  induction  hypothesis  know  that  it  is  n  —  1  times  differentiable  at  xq.  Thus,  we  can 
expand  /  as  follows  using  the  Taylor’s  theorem  with  Lagrange’s  form  for  the  remainder  term. 


/(xo  +  h) 


n— 1 


E 


r{xo) 

i\ 


h‘  + 


n! 


(73) 


for  0  <  5  <  h,  where  we  have  used  the  assumption  of  the  theorem  that  exists  in  the  neigh¬ 
borhood  around  xo,  exclusive  of  Xq. 

We  now  substitute  eq.(73)  into  eq.(12)  and  use  from  Lemma  A.l  the  fact  that  P„(h)  =  P^_i(h)  + 
Pnh^  to  get  the  following.  There  exist  constants  A  and  /iq  >  0,  and  a  polynomial  of  order  n, 
P„(x),  such  that  for  0  <  /i  <  /lo 

I  (74) 


Since  we  know  from  Lemma  A.l  that  /  is  Lipschitz  n,  from  our  induction  hypothesis  /  is  n  —  1 
times  differentiable  at  xo  and  in  fact  P„_i(h)  is  equal  to  the  first  n  Taylor  terms  of  /  about  xq. 
Thus,  eq.(74)  simplifies  to 

-  pnh^\  <  (75) 

Dividing  both  sides  of  eq.(75)  by  h",  we  get 


J'^ixo  +  S) 


n\ 


-pn\  <  Ah 


(76) 


For  any  c  >  0,  we  can  set  Ah  <  e,  whereby  for  any  S  <  h  eq.(76)  holds.  Thus,  lim5_,^o 
indeed  exists  and  is  in  fact  equal  to 


□ 
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